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ABSTRACT 

(Distribution  Limitation  Statement  No.  2) 


The  utility  (accuracy,  speed,  etc.)  of  a  hydrocode  in  solving  the  equations  of 
hydrodynamics  may  be  estimated  by  applying  the  hydrocode  test  problems  described 
in  this  report.  Given  the  estimations  of  the  utilities  of  a  pair  of  hydrocodes, 
we  may  decide  which  is  the  preferred  hydrocode  to  use.  The  hydrocode  test 
problems  described  in  this  report  have  solutions  which  are  known  exactly.  Seven 
hydrocode  test  problems  involving  shocks,  rarefactions,  and  interactions  are 
investigated  and  applied  to  a  typical  hydrocode,  AFWL*s  PUFF. 
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SECTION  I 
INTRODUCTION 


1 .  THE  HYDROCODE  UTILITY  PROBLEM 

Hydroccdes  are  computer  programs  used  to  solve  the  equations  of  hydro¬ 
dynamics.  By  “equations  of  hydrodynamics"  we  mean  the  equations  arising  from 
the  conservation  laws  and  thermodynamic  laws.  (See  Appendix  I  for  a  detailed 
discussion  of  these  equations.) 

In  particular,  we  shall  restrict  the  hydrocodes  in  this  report  to  Jbe  digital 
computer  programs  of  finite  difference  schemes.  We  shall  also  restrict  the 
discussion  in  this  report  to  one-dimensional  linear  geometry  hydrocodes.  At 
the  end  of  this  section  we  shall  give  an  example  of  a  typical  hydrocode.  The 
purpose  of  this  report  is  to  present  the  development  of  a  method  of  determining 
the  utility  (accuracy,  speed,  etc.)  of  a  hydrocode. 

Hydrocode  solutions  to  the  equations  of  hydrodynamics  usually  differ  from 
the  exact  solutions.  If  the  hydrocode  solutions  converge,  this  difference  may 
be  reduced  by  refining  the  mesh  of  the  finite  difference  scheme  used  in  the 
hydrocode.  Thus,  if  the  hydrocode  solutions  converge,  the  accuracy  of  the 
hydrocode  solutions  is  directly  related  to  the  computational  effort. 

However,  we  have  no  proof  of  convergence  of  these  schemes  except  when  the 
solution  is  assumed  smooth  between  a.  priori  known  shock  positions.*  But  the 
problems  of  interest  involve  shocks  whose  positions  are  not  a  priori  known. 
Therefore,  in  the  problems  of  interest,  we  have  no  proof  of  convergence. 

Tha  accuracy  of  the  hydrocode  solution  is  not  necessarily  directly  related 
to  the  computational  effort.  That  is,  the  mesh  (space  differences  and  time 
differences)  might  be  refined  and  yet  the  accuracy  might  not  be  improved. 
Moreover,  even  if  convergence  were  proved  for  a  scheme,  it  still  might  not  be 
of  any  practical  use  because  the  convergence  might  be  too  slow.  That  is,  the 
computational  effort  required  might  be  too  much  for  practical  purposes. 


*Lax  and  Keller,  The  Initial  and  Mixed  initial  and  Boundary  Value  Problem  for 
Hyperbolic  Systems,  Los  Alamos  Report  LAMS-1205,  T951. 
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For  the  hydrocode  to  be  useful,  the  error  should  be  small  and  the  computa¬ 
tional  effort  must  be  within  the  capacity  of  the  computer  in  relation  to  both 
memory  and  computation  time.  We  wiil  measure  the  computational  effort  by  the 
running  time  required  on  the  Air  Force  Weapons  Laboratory's  CDC  6600  computer. 

It  is  unfortunate  that  this  running  time  also  depends  somewhat  on  the  other 
programs  which  are  being  run  in  the  66u0  because  of  its  parallel  processing. 
However,  this  was  the  only  computing  time  number  available. 

When  one  uses  a  hydrocode  to  predict  some  phenomenon,  he  would  like  to 
know  how  accurate  he  can  expect  the  hydrocode's  prediction  to  be.  That  is,  he 
would  like  to  have  a  number  or  numbers  which  indicate  to  him  what  the  deviations 
between  the  hydrocode  numbers  and  the  observed  data  will  be.  One  place  for 
deviations  to  enter  is  between  the  mathematical  model  of  the  phenomenon  and 
the  phenomenon  itself.  This  report  does  not  discuss  that  problem.  (However, 
Appendix  I  does  investigate  the  mathematical  model.) 

The  problem  wa  are  trying  to  deal  with  is  the  deviation  between  the  mathe¬ 
matical  model  (represented  by  the  equations  arising  from  the  conservation  laws 
and  thermodynamic  laws)  and  the  hydrocode  model  (that  is,  the  computed  approxi¬ 
mation  to  the  solution  of  the  equations  arising  from  the  conservation  and 
thermodynamic  laws).  Therefore,  we  are  going  to  call  the  solution  to  the 
mathematical  model  the  "exact  solution."  We  therefore  want  some  number  or 
numbers  that  indicate  the  "distance"  between  the  hydrocode  solution  and  the 
exact  solution.  The  numbers  will  be  like  the  maximum  deviation  between  the 
exact  pressure  profile  and  the  hydrocode's  pressure  profile,  or  the  sum  of  the 
squares  of  the  pressure  deviations,  etc.  (See  Appendix  III  for  details  on  the 
error  numbers.)  These  numbers  will  then  give  a  measure  of  the  utility  of  a 
hydrocode  for  a  certain  type  of  problem. 

Hence,  we  developed  a  series  of  hydrocode  test  problems  which  are  represent¬ 
ative  of  the  basic  types  of  flew  encountered  in  hydrodynamic  problems  and  whose 
solutions  are  known  exactly.  By  basic  types  of  flow  we  mean  flows  including 
shocks,  rarefactions,  etc.  (For  details,  see  Section  II.)  Then  we  applied 
the  hydroc.ode  test  problems  to  a  typical  hydrocode  to  give  an  example  of  the 
procedure  that  we  have  developed  to  measure  the  utility  of  a  hydrocode. 
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2.  AN  EXAMPLE  OF  A  TYPICAL  HYDROCODE 

As  an  example  of  the  type  of  hydrocode  that  we  are  referring  to,  we  shall 
present  a  sketch  of  AFWL's  PUFF  hydrocode  which  is  described  in  detail  in 
AFWL-TR-66-48 . 

PUFF  is  a  modification  of  the  hydrocode  proposed  by  von  Neumann  and  Richt- 
myer  in  their  March  1950  Journal  of  Applied  Physics  paper.  This  paper  intro¬ 
duced  the  notion  of  artificial  viscosity.  The  equations  they  arrived  at  from 
the  conservation  laws,  thermodynamic  laws,  and  the  introduction  of  an  artificial 
viscosity  were  the  following: 

Consider  a  one-dimensional  fluid  motion.  Let  x  be 
the  Lagrangian  coordinate,  and  X=X(x,t)  be  the  Eulerian 
coordinate.  That  is,  X(x,t)  gives  the  position,  at 
time  t,  of  a  fluid  element  that  was  initially  at  posi¬ 
tion  x. 

Let  p  (x)  be  the  initial  density,  so  that  V  and  v, 
given  by 

V(x,t)  =  (l/p0(x)(3X/3x))  (1) 

and 

v(x,t)  *  3X/3t,  (2) 

are  the  specific  volume  and  fluid  velocity,  respectively. 

The  equations  of  motion,,  of  energy,  and  of  continuity 
are: 

Po(3v/3t)  =  (3/3x) (p+q) ,  (3) 

(3e/3t)  +  (p+q) (3V/3t)  *  0,  (A) 

and 

p  (3V/3t)  =  (3v/3x) .  (5) 

o 

In  these  equations,  p  =  p(x,t)  is  the  ordinary  (or 
static)  fluid  pressure  and  e  =  e(x,t)  is  the  internal 
energy  per  unit  mass.  A  connection  between  e,p,V  is 
established  by  an  equation  of  state,  which  will  be 
assumed,  for  the  purpose  of  illustration,  to  have  the 
form 

e  *  (pV)/(y-l)  (6) 

which  holds,  for  example,  in  the  case  of  a  perfect  gas. 
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Y  is  a  constant  >1. 

It  is  supposed  that  iho  dissipative  mechanism  can  be 
represented  by  the  additional  term  q  in  the  pressure, 
which  is  assumed  to  be  negligibly  small,  except  in  the 
neighborhood  of  the  shock. 


The  original  q  (artificial  viscosity)  was 


[pocAxJ  3V  3V 

V  3t  3t 


Using  (5),  it  can  also  be  written 

2 

(cAx)  3U  3U 

q  =  _ - __  _  (8) 

where  c  is  a  dimensionless  constant  near  1  and  Ax  is  the  interval  length  used 
in  the  numerical  integration  of  the  hydrodynamic  equations.  PUFF's  q  is  a 
modification  of  this. 


Now  we  will  describe  PUFF’s  difference  scheme  to  solve  the  foregoing  equa¬ 


tions. 


Let  the  points  of  a  rectangular  network  with  spacings  Ax  and  At  be  denoted 
by  x  ,  tn  (£  *  0,1,2,  ...,  L;  n  =  0,1,2,  ...  ).  We  shall  also  have  occasion 

X* 

to  deal  with  intermediate  points,  having  coordinates 


£+1/2 


=  1/2(x£+1+x£),  tn+1/2  =  l/2(tn+1+tn). 


To  facilitate  the  writing,  we  introduce  abbreviations  such  as 


£+1/2  \  £+1/2 


vk+i/2-t”) etc- 


The  difference  equations  with  which  FUFF  approximates  the  differential 
equations  (1)  through  (8)  are  the  following: 

n+1/2  n-1/2  n  n-1/2  pn  n-1/2 

U£  ~  U£  =  _  *£+1/2  q£+l/2  *£+1/2  q£-l/2  ... 
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is  PUFF{s  difference  approximation  to 

au^.MP+s ) 

po  9t  5x 


(10) 


Yn+1  yn 

X-lh.  m  tft+1/2 

At  5, 


is  PUFF's  difference  approximation  to 


3X 
9 1 


tt 


V 


n+1 

p£-l/2 


PoAx 


Yn+1 

i 


.n+1 

"£"1 


(ID 


is  PUFF's  difference  approximation  to 


po  9X  V  9X 

p  ”  9x  °r  V  9x 

o 


because 


P 


Now  let 


AU  =  yn-rl/2  yn+1/2 


Then  PUFF's  q  is 


n+1/  2 
q £-1/ 2 


(12) 
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where  CQ  =  1.8,  =  0.25,  Cs  =  isothermal  sound  speed  at  time  tn and 

position  £-1/2.  That  is, 


s  dp 


e  const. 


Then  PUFF  simultaneously  solves 


n+i  _  n  n+1  n+1/2  n 

n  t-1/2  £-1/2  *£-1/2  *  q£-l/2  *£-1/2 

At  2 


n-1/2 

q£-l/2 


Pn^  s  P 
* £-1/2  P 


/m-1  n+1  \ 

y  £-1/2*  p£-l/2/ 


(the  equation  of  state). 

The  difference  equation  is  the  result  of  differencing 


n  ,  /n,  n  1  3D 

0  at  + 


which  results  from 


n  x  9V  .  3V  1  3U 

0  *  aT+  <pT’>  aland 


PUFF's  method  of  solutiovi  is  this:  Suppose  that  all  quantities  are  known 

for  superscript  n  or  n-1/2  (this  is  referred  to  ac  being  at  cycle  n );  compute 

for  each  £  from  (9),  then  compute  for  each  £  from  (1);  next  compute 

n+1/2 

p£-l/2  ^or  eac^  *  ^rom  then  compute  q^_^/2  ^or  sach  £  from  (12);  next 

•  n*f*l  * 

compute  ££_-.y2’  **£-1/2  ^or  eac^  ^  s^jnu^taneous^y  solving  (13)  and  (14).  At 

this  point,  all  variables  have  been  advanced  up  to  cycle  n+1.  Next,  PUFF  does 

its  time-step  computation  to  compute  the  next  At.  The  time-step  computation 

is  based  on  a  modification  of  the  criterion  developed  by  Courant,  Friedrichs, 
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and  Levy  in  their  1928  MATH.  ANNALEN  paper.  The  physical  interpretation  of 

this  criterion  is  that  the  time  step  should  be  restricted  so  that  a  sound 

signal  cannot  travel  across  more  than  one  zone  in  a  time  step.  That  is, 

A*-1*  £  f°r  all  £>  where  is  the  sound  speed  in  the  £th  zone  at  cycle 

n,  At  is  the  nth  time  step,  and  AX1^  is  the  width  of  the  £th  zone  at  cycle  n. 

Thus,  Atn  is  usually  defined  as 

Yn 

At  s  0  min  -- 

where 

0  <  0  <  1 

0  is  called  the  CFL  number.  PUFF  uses  a  modification  of  this  (see  AFWL-TR- 
67-48).  For  more  details  about  the  CFL  criterion,  see  Appendix  II 
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SECTION  II 

DESCRIPTIONS,  DERIVATIONS,  AND  DETAILS 

To  make  t:his  section  easier  to  read,  we  include  a  discussion  of  the  format 
used  in  presenting  the  hydrocode  test  problems.  Seven  hydrocode  test  problems 
are  considered.  The  problems  discussed  in  this  report  are  labeled  SCTP  for 
Slab  Code  Teat  Problems.  The  "Slab"  refers  to  the  geometry,  i.e=»,  one- 
dimensional  linear. 

The  labeling  is  further  modified  to  indicate  which  problem,  i.e.,  SCTP-I, 
SGTP-II,  ...»  SCTP-VII.  The  consideration  of  each  problem  is  arranged  in  the 
following  way: 

n.  HYDROCODE  TEST  PROBLEM  SCTP-n 

a.  Description  of  Problem 

Here  we  give  a  graphical  and  verbal  description  of  the  problem. 

b.  Derivation  of  Solution 

Under  this  heading  the  exact  mathematical  solution  is  derived. 

c.  Application  as  a  Test  Problem 

Now  we  get  down  to  the  numerous  details  and  difficulties  of 
applying  the  problem  to  a  hydrocode,  that  is,  how  one  inputs 
the  initial  and  boundary  values  for  various  codes.  The  first 
division  in  variety  of  codes  is  whether  the  code  is  Eulerian 
or  Lagrangian  in  its  formulation  of  the  hydrodynamics  equa¬ 
tions.  Therefore,  we  introduce  the  following  subdivisions. 

(1)  Eulerian  Input 

(a)  Initial  values 

(b)  Boundary  values 

(2)  Lagrangian  Input 

(a)  Initial  values 

(b)  Boundary  values 
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In  the  foregoing  subdivisions  we  describe  how  the  hydrocode  test  problems 
may  be  introduced  into  either  an  Eulerian  or  a  Lagrangian  hydrocode.  That  is, 
we  give  the  specific  processes  for  inputting  the  initial  values  and  the 
boundary  values. 

Up  to  this  point  we  have  not  specified  any  numbers  for  the  initial 
pressures,  densities,  etc.  These  numbers  are  introduced  in  the  following 
subdivision. 

(3)  Numerical  Values  for  SCTP-n 

Here  we  give  the  initial  and  boundary  value  numbers  which 

are  used  for  both  Eulerian  and  Lagrangian  hydrocodes. 

Also,  we  give  the  resulting  numbers  which  occur  In  the 

solution. 

Next,  we  discuss  what  is  to  be  expected  from  these  problems  when  they  are 
applied  as  hydrocode  test  problems.  This  is  done  under  the  heading 

(4)  Comments  on  the  Computer  Solution 

Under  this  heading  we  have  two  subdivisions: 

(a)  General  comments 

In  this  paragraph,  we  give  a  discussion  of  what  can 
be  said  about  the  hydrocode  test  problems  behavior 
in  general.  That  is,  how  long  the  solution  is  as 
described  without  altering  the  boundary  values, 
what  mistakes  a  Lagrangian  hydrocode  will  probably 
make,  what  errors  an  Eulerian  hydrocode  will  probably 
commit,  etc. 

Lastly,  as  a  final  illustrative  help  to  the  reader  who  wishes  to  apply  the 
hydrocode  te3t  problems,  we  present  an  application  of  the  test  problems.  We 
apply  the  seven  hydrocode  test  problems  to  the  typical  Lagrangian  hydrocode 
PUFF  (described  in  the  introduction) .  This  is  done  in  the  last  heading: 

(b)  PUFF  comments 

Here  we  present  the  error  graphs  and  error  tables 
for  PUFF  to  illustrate  what  we  believe  to  be  a 
reasonable  way  of  describing  the  "distance" 
between  the  PUFF  solutions  and  the  exact  solutions. 
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We  did  not  illustrate  the  Eulerian  input  by  applying  the  hydrocode  test 
problems  tc  an  Eulerian  cede. 

1.  HYDROCODE  TEST  PROBLEM  SCTP-I 

a.  Description  of  Problem 

This  is  the  steady  profile  solution  of  a  constant  velocity  piston 
compressing  the  fluid  ahead  of  it.  By  steady  profile  solution  it  is  meant 
that  the  fluid  parameters  from  the  piston  face  to  the  shock  front  are  constant, 
and  also  to  the  right  of  the  shock  front  the  fluid  parameters  are  constant. 
Therefore,  the  velocity  of  the  piston  and  the  velocity  of  the  fluid  between  the 
piston  and  the  shock  front  are  the  same. 

Figure  1  attempts  to  explain  the  problem  graphically. 


Figure  la.  Pipe  Plot  for  SCTP-I 
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v^.  5  fluid  velocity  to  the  right  of  the  shock 

v  =  shock  velocity 
s 

b .  Derivation  of  Solucion 

(1)  Conservation  of  Hass  or  the  Equation  of  Continuity 

Looking  from  the  right  side  of  the  shock,  wave,  the  mass  engulfed 
in  time  At  by  the  shock  wave  is 


(v  -v  ^  p  AAt 
s  cj  r 

Looking  from  the  left  side  of  the  shock  wave,  the  mass  cowing  through  the 
shock  wave  in  time  At  is 


(vvi)  "i 

Let  m  be  the  mass'  per  second  per  unit  area  passing  through  the  shock.  Then  in 
time  At,  the  mass  mAAt  will  pass  through  the  shock-  Hence, 


(1:00 


So, 


v£  *  v.  -  m/p£  and  vr  »  vg  -  a/pr 


Tharefores 


(2)  Newton’s  Second  Law  or  the  Equation  of  Hotien  or  Conservation  of 
Momentum  Applied  to  Hass  Passing  through  the  Shock  Wave 

Suppose  the  mass  mAAt  passes  through  the  shock  wave  ia  thee  At. 
From  F  *  d/dt  (Mv)  comes 


lin 

At-*C 


(cAAt)v£  -  (oAAt)vr 
At 
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So, 


F  =  niA 


(v!-°r) 


Also, 


F  = 


(prFr)A 


Therefore, 


h  -  Pr  *  m(vrvr) 


(2tN2L) 


(3)  Conservation  of  Energy  or  Energy  Balance  Equation 
(a)  Hark  done  on  a  mass  passing  through  shock  wave 
Recall  that 


V  *  J 


F  -  dS 


The  work  dene  on  the  mass  mAAt  entering  the  shock  wave  is  ^v^Atj- 

The  work  done  on  the  mass  mAAt  exiting  the  shock  wave  is  ^v£At^. 

Therefore,  the  work  dona  on  the  mass  mAAt  as  it  passes  through  the  shock  wave 


xs 


(Vs '  Vi)“l 


AH  =  the  work  done  on  the  mass  mAAt  per  area  A  per  time  At. 


AH  =  ?  ✓„  -  P  v 
i  l  r  r 


(b)  Kinetic  energy  change  in  a  mass  passing  through  the  shock 


wave 


The  change  in  kinetic  energy  of  the  mass  mAAt  as  it  passes 
through  the  shock  wave  is  1/2  mAAt 

AK  =  tne  rhange  in  kinetic  energy  of  the  mass  mAAt  per  area  A  per  time  At. 
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AK  =  1/2  m^v|-v2  j 


(c)  Internal  energy  change  of  a  mass  passing  through  the  shock 

wave 

The  specific  internal  energy  e  is  a  function  of  the  pressure 
and  specific  volume  e  =  e(P,V).  So  the  internal  energy  of  a  gas  with  mass  mAAt 
is  (mAAt)  e  (P,V). 

AI  =  the  increase  in  the  internal  energy  of  the  mass  mAAt  per  area  A  per  time  At 
as  it  passes  through  the  shock  wave  in  time  At. 


AI  -  m 


(•  e  M) 


(d)  Energy  accounting 

The  increase  in  the  mass  mAAt's  internal  energy  plus  the 
increase  in  its  kinetic  energy  is  equal  to  the  work  done  on  it.  Symbolically, 
AAt  (AIvAK)  =  (AW)  AAt  or  AI  +  AK  *  AW  or 

"  (e  (Vt)  -  e  (VVr))  +  I  (’H ) *  V*  -  Vr  (3:C® 


(4)  Equation  of  State 

The  equation  of  state  used  is  that  of  the  perfect  or  ideal  gas 
model  of  standard  air. 


(y-1)  e  =  PV  =  RT  (4: EOS) 

y  =  the  ratio  of  specific  heats.  It  will  be  taken  as  1.4  which  is 
approximately  the  value  for  standard  air 

e  =  specific  internal  energy 

P  =  pressure 

V  2  specific  volume 

R  =  gas  content 

T  =  temperature 
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(5)  The  Rankir.e-Hugoniot  Relation 


By  (2:N2L) 


By  (1:CM)  ' 


By  (4: EOS) 


(VPr)  (vvr)  "  n  (vvr)  (Vvr) 

f  H^r)  ‘  (Vr)^) 

n  =  (vv r)  '  (Vr~V£ ) 
e  (P,V)  *  PV/(y-l) 


Substituting  these  into  (3:CE)  yields 


£r  (v*.) +  (prpr)  “ V*  -  V. 


0  *  e  ~  e 
%  r 


n+P  ,  , 

(v\; 


(5:RH) 


(6)  The  Shock  Speed  Relations 
Recall  that  by  (1:CM) 


(VVt)  ‘  (Vvr)  p, 


By  (2:N2L) 


=  m  (vvr) 


By  (3:CE) 


-2—  fp  V  — p  V  \  +  -Jr  \  =  P„v„  -  P  v 

y-1  \  i  l  r  ij.  l  \  i  r)  l  i  r  r 
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Eliminating  m,  p^,  V^,  and  P  ,  these  equations  yield 


1 

Y-l 


P  + 
r 


(v  -v  )  p  (v.-v  — — —  V  -  P  V 

\  s  r/  "r  \  i  r/j  vg-vr  r  r  r 


This  reduces  to 


K+-I 

[v  -V  ) 

Ur/ 

pr  (7lfvr) 

'  ?rvr 

/V  "V  \  p 

lS  X)  1 

The  solution  to  this  is 


(6:SS(v)) 


The  plus  sign  is  taken  on  the  radical  so  that  vg  >  0. 
From  (1:CM)  and  (2:N2L) 


P»  '  'r  ’  (V\)  (V\) 


P  ~P 
%  r 


=  v  v.  - 
s  Si 


V  V 

s  r 


-  v  v6 
r  a 


+ 


Therefore 
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Substituting  this  into  the  v  quadratic  yields 

s 


V  V 
s  r 


=  0 


or 


o 

v 

s 


-  2v  v 
r  s 


+  v2 


o  _p 

-  J±l2A_£ 


-  O  = 
r 


So, 


v 

s 


p  _p 

1+Y  *  rr 


+  C2 


So 


» 


v 

s 


V 

r 


By  (5:RH) 


Substituting  this  into  (?:SS(P))  yields 


(7:SS(P)) 


v 

s 


2yP  V2 
r  r 


(y+1)  v£  -  (v-i)  vr 


(8:SS(V)) 
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Therefore,  for  y  =  1.4,  p./p  <  6  Is  required. 

J6  r 

(7)  Ad litional  relations 

Suppose  the  right  side  values  are  specified  and 

(a)  is  specified 

By  (6:SS(v)),  v  is  determined, 
s 

By  (1:00  and  (2:N2L) 

’  Pr  +  (V'r)  pr  (V'r) 

By  (1:00 
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(c)  v£  is  specified. 

Recall  (8:SS(V))  is 


+  /  2^PrVr 

vs  "  vr  J  (y+1)  V£  -  (y-1)  Vr 


In  order  for  the  radical  to  be  defined,  it  is  required  that 


V  <  V 

i  r  y+1 


If  this  restriction  is  satisfied,  then  v  is  determined. 

s 


By  (1:CH) 


V£  Vs 


-  V  (v  -V  ^  p 
l  \  s  r/  *r 


By  (2:N2L) 


p  =  p 

P£  *r 


+  (V  -v  }  p  |v  -V  } 
\  s  r/  r  \  f,  r/ 


Lcatlon  as  a  Test  Problem 


(1)  Eulerian  Input 

(a)  Initial  values 


To  the  right  of  X  ,  the  initial  shock  position,  define  the 
s 

velocity  to  be  v^  *  0,  the  pressure  to  be  >  0,  the  density  to  be  pf  >  0. 
From  this,  c£  =  y determined. 

Let  Vp  =  MC^  for  some  M  >  0.  Now  the  shock  velocity  is 
determined  by  (6:SS(v)), 


vs  =  CrM  Vti  l1  + J  1+  TFITh 
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For  0  £  X  <  Xg,  the  velocity  is  defined  by  =  v  ,  the  pressure  is  defined  by 

?  =  P  +  v  v  p  .  (2:N2L)  and  (1:CM)  defines  the  density  p„  =  v  p  / fv  -v 
*  r  s  P  r  sKr'  \  s  p/ 

(b)  Boundary  values 

At  X  =  0,  hold  the  velocity  at  v^,  the  pressure  at  P^.  and 
the  density  at  p^.  Notice  that  for  an  Eulerian  code,  this  is  as  if  the  piston 
were  always  to  the  left  of  0  and  the  computer  solution  were  exact  from  the 
piston  to  0. 

(2)  Lagrangian  Input 

(a)  Initial,  values 

The  Lagrangian  input  is  the  same  as  the  Eulerian  input. 

(b)  Boundary  values 

The  Lagrangian  input  is  the  same  as  the  Eulerian  input. 
(Notice  that  the  interpretation  of  the  boundary  value  is  much  more  satisfactory 
now  because  the  piston  is  always  at  x  =  0.) 

(3)  Numerical  values  for  SCTP-I 
(a)  SCTP-I- A 

M  =  1,  AX  =  1  meter,  X  =  50  meters 
P^  =  104  dynes/cm2 
Pr  *  10"°  gm/cm3 


Right  boundary  at  300  meters.  These  values  yield 
v^  ^  1.18  x  105  cm/sec 
vs  ~  2.08  x  105  cm/sec 
P^  ^  3.47  x  104  dynes/cm2 
4.34  x  105  cm3/gm 

A  reasonable  output  recipe  is  to  run  the  problem  to  0.1  second  with  prints  at 
0.01  second  intervals.  The  300r.h  zone  should  remain  inactive  so  total  energy 
sums  will  be  taken  out  to  there.  With  the  CFL  number  set  to  1,  the  first  At 
will  be  about  6.9  x  10“4  sec.  After  the  first  time  step,  the  time  steps  should 
all  be  about  3  x  10_4  sec.  Therefore,  it  should  take  about  335  cycles  for  a 
CFL  number  of  1. 


20 


'-V'i'V* 


AFWL- TR-6  7-127 


<b)  SCTP-I-B 

Same  as  A  except  M  =  100.  This  yields 

~  1.18  x  IQ7  cm/sec 

—  1.68  x  108  dynes/cm2 

V  cz  1.67  x  105  cm3/gm 

~  6.26  x  108  cm/sec 

v  a  1.42  x  107  cm/sec 

s 

A  reasonable  output  recipe  is  to  run  the  problem  to  1C“3  sec  with  prints  at 
10"4  sec.  With  a  CFL  number  of  1,  the  first  time  step  will  be  m;  1.6  x  10~5  sec. 
After  the  first  time  step,  the  time  steps  should  run  about  2.7  x  10~6  sec. 

Thus,  it  should  take  about  375  cycles  to  run  to  10“  3  sec  if  the  CFL  number  is  1. 
.As  in  SCTP-I-A,  the  3G0th  zone  should  remain  inactive  so  total  energ/  sums  will 
be  taken  from  zone  1  to  zone  300. 

(4)  Comments  on  the  Computer  Solution 

(a)  General  comments 

If  the  hydrocode  is  solving  this  problem  correctly,  the 

velocity  of  the  shock  will  be  about  vg .  The  quantities  to  the  left  of  the 

shock  should  remain  at  v^,  P^,  p^,  and  the  quantities  to  the  right  should  be 

v  ,  P  ,  p  . 
r»  r>  '-r 

The  shock  front  should  remain  sharp  and  not  smeared  over  too 
many  zones.  The  specific  internal  energy  on  the  left  should  be  P^V^/(y-l);  on 
the  right,  it  should  be  PrVr/(Y-l).  The  specific  kinetic  energy  on  the  left 
should  be  1/2  v2;  on  the  right,  it  should  be  1/2  v2. 

The  specific  total  energy  should  be  1/2  v2  +  P^V^/Cy-^l)  on 
the  left  and  PrVr/(y-l)  on  the  right.  To  check  conservation  of  total  energy/ 
sum  ^1/2  v^  +  P^Vj/(y-])j  Mj.,  where  M^.  is  the  mass  of  the  Ith  zone,  from  the 
piston  face  to  some  fixed  point  on  the  right  which  remains  undisturbed  through¬ 
out  the  problem  (zone  300  for  SCTP-I-A  or  SCTP-I-B).  This  value  should 
increase  by  P^v^t,  where  t  is  time,  for  3  Lagrangian  code.  For  a  Eulerian 


ft 


code,  take  the  above  sum  over  I  and  subtract  t  v^l/2  v2  +  P^V^/ (y-1) j  ?} 
because  of  the  interpretation  mentioned  in  (l)(b).  After  this  is  subtracted 
from  the  total  energy  sum,  the  remainder  should  increase  by  P^v^t.  The  reason 
for  the  time  step  sequence  as  mentioned  in  (3) (a)  is  that  if  the  CFL  number  is 
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1,  then  for  the  first  time  step  At^  =  -AX/C^  ~  6.9  x  10”'1;  but  after  one  of  the 
zones  on  the  right  has  been  compressed,  it  should  be  about  43.4  cm  long;  then 

at,  -  10-*  sec 

C1 

The  shock  speed  is  moving  °.t  about  2  x  105  cm/s ec,  so  in  the  first  time  step 
it  moves  about  140  cm,  which  will  engulf  the  first  zone  to  the  right.  After 
the  first  time  step,  the  first  zone  on  the  right  should  be  compressed  to  about 
43.4  cm. 

Care  should  be  taken  in  determining  the  correct  time  ant; 
position  of  all  quantities  produced  by  the  code  being  tested.  For  example,  in 
PUFF,  the  velocities  are  one  half  time  step  behind  in  time  and  the  Jth  values 
of  the  pressures,  densities,  and  internal  energies  are  halfway  between  X(J-l) 
and  X(J).  Errors  on  the  order  of  5  or  10  percent  can  arise  if  these  variables 
are  not  plotted  at  rhe  right  time  and  place.  The  energy  partition  is  cf 
interest  because  it  indicates  whether  or  not  the  energy  dissipation  rate 
associated  with  entropy  increase  across  a  shock  is  correct.  In  particular*  if 
the  code  uses  an  artificial  viscosity  as  a  dissipation  mechanism,  the  energy 
partition  can  indicate  whether  It  introduces  too  much  or  too  little  dissipation. 

(b)  PUFF  Comments 
1.  SCTP-I-A 

ACCURACY:  Overall  accuracy  was  on  the  order  of  1  percent 
except  for  the  region  of  the  initial  discontinuity  zone  and  the  zones  where  the 
shock  transition  occurs.  That  is,  P,  V,  v,  e,  and  the  shock  speed  were  all 
within  about  1  percent  except  at  zone  52  (zone  directly  in  front  of  the  initial 
discontinuity)  where  e  and  V  were  about  13  percent  high  and  in  the  shock  transi¬ 
tion  zones  where  the  P,  V,  v,  e  take  about  six  zones  to  change  from  within  1 
percent  of  Lite  left  ^tete  to  within  .1  percent  of  the  right  state.  To  change 
from  10  percent  of  the  left  state  to  10  percent  of  the  right  state  takes  about 
three  zones.  I  believe  that  the  e,  V  peak  in  zone  52  was  caused  by  the  conver¬ 
sion  pf  kinetic  to  internal  energy  by  the  artificial  viscosity  which  was 
intensified  here  by  the  sharp  jump  at  the  initial  shock  wave  position.  After 
the  shock  had  been  rounded  off  a  bit,  the  artificial  viscosity  term  reduced 
and  remained  uniform  for  the  balance  of  the  problem.  This  increase  in  internal 
energy  combined  with  constant  pressure  produced  a  higher  specific  volume  because 
the  specific  volume  is  proportional  to  the  internal  energy/pressure  quotient. 
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Beh’ud  the  shock,  there  is  a  little  rarefaction  dip  in 
density  and  pressure  which  is  traveling  to  the  left.  At  early  times  this  dip 
exceeds  1  percent  but  does  n>t  reach  2  percent,  and  at  later  times  it  fades 
away. 

The  above  statements  indicate  that  the  Rankin e-Hugoniot 
relations  are  being  satisfied  asymptotically.  PUFF’s  total  energy,  total 
internal  energy,  and  total  kinetic  energy  errors  (deviations  froa;  the  exact 
solutions)  are  small  as  indicated  in  Table  I,  thus  indicating  approximate 
conservation  of  energy  and  correct  energy  partition.  In  particular,  the  Q 
(artificial  viscosity)  conversion  of  kinetic  energy  to  internal  energy  was 
operating  at  about  the  correct  rate. 

TIME:  PUFF  took  1033  cycles,  giving  it  an  "effective 
CFL  number"  of  335/1033  «  0.32.  PUFF's  rezont  capability  was  not  utilized. 

FUFF  took  72  seconds  CDC  66GG  central  processor  time  to  run  1033  cycles  on 
this  problem. 

2.  SCTP-I-B 

ACCURACY:  Again,  the  only  distortions  were  associated 
with  the  initial  discontinuity  and  the  later  shock  posivions.  The  internal 
energy  peak  in  zone  52  was  at  about  12.7  x  1013  ergs/gm  wrxle  it  should  have 
been  about  7.00  x  1013  ergs/gn.  The  pressure  was  correct  it  zone  52  so  there 
was  a  corresponding  peak  in  the  specific  volume  at  zone  52.  Vhe  rarefaction 
dip  in  the  pressure  and  density  which  is  seen  at  early  tires  is  less  than  3 
percent.  By  cycle  1495,  this  rarefaction  dip  is  down  to  1  percenv  or  less. 

The  shock  transition  is  again  about  six  zones  for  1  percent  of  the  left  to  1 
percent  of  the  right  and  about  three  zones  for  10  percent  to  10  percent. 

Except  for  the  above  mentioned  distortions,  the  PUFF  solution  to  SCTP-i-B  has 
at  most  about  I  percent  error  for  all  other  zones.  Again,  the  sum  total  energy, 
sum  internal  energy,  and  sum  kinetic  energies  are  close  as  indicated  by  Table  I. 

TIME:  PUFF  took  1463  cycles,  giving  an  effective  CFL 
number  of  375/1463  »0.26.  PUFF  took  75  seconds  CDC  6600  central  processor 
time  to  run  1463  cycles  on  this  problem. 

2*  Discussion  of  the  initial  discontinuity 

In  inputting  a  discontinuous  velocity  profile  into  PUFF, 
there  are  several  alternatives  in  selecting  v  (the  value  of  the  velocity  at 
the  discontinuity).  Graphics? ly,  the  exact  picture  is  this: 
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v 


where  m  Is  the  mass  in  the  zone  to  the  left  of  the  discontinuity  and  m  is  the 
£  r 

mass  of  the  zone  to  the  right  of  the  discontinuity.  X  is  the  position  of  the 

s 

discontinuity,  Xs  l  is  the  adjacent  zone  boundary  to  the  left,  and  is  the 

zone  boundary  to  the  right. 

Since  PUFF  uses  a  linear  interpolation  on  its  velocities, 
the  PUFF  profile  will  look  like  this  when  a  value  v^  is  given: 


v 
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Same  alternatives  for  v  i 

m 

(1)  Set  v  to  either  v  or  v  . 

*-*  £  x 

(2)  Minimize  the  maximum  deviation  of  PUFT’s  velocity 
profile  from  the  exact  velocity  profile.  This  leads  to 


v£+v r 


yields 


(3)  Conserve  the  momentum  of  the  exact  profile.  This- 


m  v  +  m  v 

v  =  £  r  r 

ra  mn-+m 

?.  r 


(4)  Conserve  the  energy  of  the  exact  profile.  This  yields 


/WVr'j 

\  V"r  / 


VftV 
l  V* T  , 


/V'Hj  \ 

V  Vr  / 


which, iwhen  v_  =  0,  becomes 


\  /  °i+mr 


Notice  that  if  vr  -  0,  the  vm  of  (4)  is  bigger  than  the  v^  of  (3)  which  in  turn  is 
bigger  tnan  the  of  (2)  when  <  m^.  So  we  have,  when  v^  =  0  and  m£  <  m  , 

V_  <  V(2)  <  v(3)  <  V(4>  <  V 


Now  the  internal  energy  spike  increases  as  one  increases  the  upper  index  of 

j  (1)  y  *  • 

Vm  \~e*‘  vn  =  vr  and  v^  **  the  rarefaction  dip  and  the  associated 

hump  leading  it  decrease.  In  all  subsequent  shock  inputs  in  this  report,  we 


set  v  =  v. 
ra  l 
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Table  I -A* 

- 

PUFF  ERROR  TABLE  FOR 

SCTP-I-A  WITH  v  =  v 
m  2 

Problem  time  *» 

Computer  time  = 

0.1  sec 

72  sec 

PUFF  cycle  =  1033 

Number  of  active  zones 

=  267 

Sum  Sqr.  Frr. 

Max.  Err. 

Position  of 
Max.  Err. 

X 

Pressure 

0.0251 

+0.339 

Velocity 

0.0423 

+0.583 

s 

X 

s 

X 

s 

Density 

0.0229 

+0.292 

Sum  Int.  Energy 

Sum  Kin.  Enerpv 

Sum  Tot.  Energy 

EXACT 

1.324  x  10s 

2.270  x  108 

1.551  x  10? 

PUFF 

1.324  x  109 

2.264  x  108 

1.550  x  109 

Table  I-B 

rUFF  ERROR  TABLE  FOR  SCTP-I-B  WITH  -  =  v 

'm  l 

Problem  time  « 

3  1.0  x  10  8  sec 

PUFF  cycle  =  1463 

Computer  time 

*75  sec 

Number  of  active  zones 

=  198 

Sum  Sqr.  Err. 

Max.  Err, 

Position  of 
Max .  Err . 

Pressure 

0.0486 

+0.644 

X 

s 

X 

g 

Velocity 

0.0764 

+0.850 

Density 

0.0568 

+0.577 

X 

s 

Sum  Int.  Energy 

Sum  Kin.  Enerev 

Sum  Tot.  Energy 

6.188  x  1012 

EXACT 

3.095  x  .10 12 

3.093  x  1012 

PUFF 

3.098  x  l^12 

3.037  x  1012 

6.135  x  1012 

*See  Appendix  III  for  a  discussion  of  the  Error  Tables. 
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Equation  of  state: 


(y-1)  e  -  PV 


Since  the  process  is  homentropic. 


and 


So, 


dP 

dp 


lX 


dP 

dp  PX 


Therefore,  the  equations  form  the  homogeneous  hyperbolic  system: 


The  method  of  characteristics  will  be  used  to  derive  the  solution. 

(2)  The  Method  of  Characteristics 
(a)  Canonical  form 

For  convenience,  introduce  a  sue  that  do  =  Cdp/p.  Fo 

sy 


P  *  P 


GJ 


o  =  2C/(y-l)  +  const.  So,  define  o  as  2C/(y~1)*  Then 
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Let  v  =  v  -f  o,  to  =  v  -  0.  Then 


Notice  that  the  matrix  is  now  In  characteristic  form.  That  is,  its  character¬ 
istic  values  are  on  the  diagonal  and  zeros  are  elsewhere. 

(b)  Characteristic  curves 


Consider  the  curves  in  the  X,t  plane  which  are  solutions  of 

the  ordinary  differntial  problems: 

v+/ +\  v+  ,  dX  _ 

X  }t  }  =  X  and  -rr—  »  v  +  C 
\ o/  o  dt 

or 

X  (t  )  *  X  and  ™—  *  v  -  C 
V  o  /  a  dt 

Notice  that 


d_ 

dt 


x  dt 


+  v. 


0 


and  likewise. 


■JjT  tu^X  (t),tj  «  0 

v  and  to  are  called  Rlemann  invariants  because  v  is  constant  on  the  curve 

(x+(t),tj  which  is  called  a  characteristic  curve  and  likewise,  m  is  constant 
'  '  /  „  \ 
on  the  characteristic  curve  (X  (t),t^.  Therefore, 

v(X+(t)>t)  .  v(x+,tj) 

and 

<o(x"(t),t)  -  u(x“,t") 

Hence,  if  v^Xo,toj,  a^XQ,toj  are  known  for  some  (xQ,to)  the  values  of  v  +  o 
are  known  at  all  points  along  the  curve  ^X— (t),tj.  The  initial  and/or  boundary 
data  for  the  original  partial  differential  problem  car.  be  used  to  determine 
values  of  v  and  o  at  some  fxo,tV  That  is,  the  initial  data  specifies 
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v(xQ,t^j  and  o^XQ,t for  tQ  =  0.  We  repeat  that  the  curves  ^X*(t),tj  and 
|x  (t),t^  are  called  characteristic  curves.  The  family  of  curves  in  the  (X,t) 
plane  with  slope  v  +  C  is  designated  C+  characteristics.  Likewise,  the  family 
of  curves  in  the  (X,t)  plane  with  slope  v  -  C  is  designated  C_  characteristics. 
Notice  that  on  a  C+  characteristic  dv/dp  *  dv/dp  +  do/dp  =  0  and  do/dp  =  C/p 
so  dv/dp  =  -  C/p.  With  this  relation,  we  may  define  a  mapping  from  C+  charac¬ 
teristics  into  curves  in  the  (v,p)  plane.  Let  the  curve  ^X+(t),t^  with  initial 
values  ^XQ,to^  be  mapped  into  the  curve  ^v(p),p^  with  initial  values  (V0»P0^ 
where  vq  =  v^0*t0^»  P0  =  P  (XQt  tQ^  and  dv/dp  =  -  C/p.  This  family  of  curves 
in  the  (v,p)  plane  is  designated  characteristics  (because  they  are  the 
images  of  the  C+  characteristics).  Likewise,  the  family  of  curves  in  the 
(v,p)  plane  with  slope  C/p  is  designated  characteristics  (because  they  are 
the  images  of  the  C  characteristics). 


The  method  of  characteristics  will  now  be  applied  to  SCTP-II. 
Remember  that  SCTP-II  can  be  interpreted  as  a  piston  withdrawing  from  a  pipe 
containing  a  gas  initially  at  rest.  Let  Xp(t)  be  the  position  of  the  piston 
at  time  t,Xp(0)  «  0.  Let  vp  be  the  velocity  of  the  piston,  vp  <  0.  Now  the 
information  that  the  piston  is  being  removed  travels  back  into  the  undisturbed 
gas  with  the  speed  of  sound  of  the  resting  gas,  C .  So,  in  the  X,t  plane  a 
wave  plot  looks  like  figure  2b.  First  C^Xp(t),t^,  the  sound  speed  at  the  face 
of  the  piston,  will  be  determined  by  the  method  of  characteristics.  Notice 
(in  figure  2b)  that  the  C_  characteristics  have  ^Xo,tQj  in  the  undisturbed 
region  (therefore,  v^XQ,to^and  a(xQ,to^  are  known)  and  the  C_  characteristics 
intersect  the  piston  path  line.  The  C_  characteristics  will  intersect  the 
piston  path  line  because  a  C_  characteristic  will  start  out  in  the  undisturbed 
region  with  slope  -Cr;  than  where  the  C_  characteristic  enters  the  rarefaction 
wave,  its  slope  v  -  C  becomes  more  negative  than  the  fluid  velocity  as  long  as 
C  >  0.  C  will  be  positive  if  the  piston  is  not  withdrawn  too  fast  for  .the 
gas  to  follow  (that  is,  if  no  vacuum  is  formed  between  the  gas  and  the  piston) . 
When  no  vacuum  is  formed,  the  gas  will  be  in  contact  with  the  piston  face  and 
there  it  will  have  the  velocity  vp  of  the  piston.  Therefore,  the  C_  charac¬ 
teristics  intersect  the  piston  path.  Hence,  if  ^XQ,r 


,o) 


is  in  the  undisturbed 


region  (e.g.,  if  tQ  «  0  and  X^  >  0),  v^XQ,to^  =  0,  o^XQ,to^  ■“  2Cr/(Y-l),  and 
v(xo,to)  -  o(XQ,tJ  =  v(xp(t),t)  -  o(xp(t),t).  So,  -2Cr/(y-l)  *  v?  -2c(xp(t) ,  t)/ 
(y-1).  Therefore,  C  Xp(t),t  »  y-1/2  vp  +  Cr. 
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Figure  2b.  Eu.lerian  Wave  Plot  for  SCTP-II 
(c)  Simple  waves 

Note  that  v  -  a  is  the  same  constant,  -o  ,  for  the  entire 

r* 

rarefaction  region.  Such  a  region  is  called  a  C_  simple  wave.  Also,  if  v  +  c 

is  the.  same  constant  in  a  region,  then  the  region  is  called  a  C+  simple  wave. 

Since  v  -  o  =  -o^  in  this  region,  v  s  u  -  so  v  depends  only  on  o.  But  o 

depends  only  on  C  which  depends  only  on  p.  Therefore,  P,  v,  a,  C  depend  only 

on  p.  Notice  that  if  o  -  c  is  substituted  for  v  in 

r 


then  the  one  equation  + 


Hr*:) 


=  0  results.  Since  C  =  y-1/2  o,  this 


equation  is  of  the  form  oc  +  (ao+b)  o^  =  0  where  a  and  b  are  constants.  The 
general  solution  to  equations  of  this  type  may  be  expressed  implicitly  as 
o(X,t)  -  F^X  -  (ao(X,t)  +  b)t). 
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So, 


?(x  -  (a-or+c)t)  =  f(x  -  (v+C) t^  =  F^X  -  (^  °'°o)r) 


o  (X,t)  =  F^ 

F  can  be  determined  by  the  condition  at  the  piston  face.  For  t  <  0, 


°(Xp(t)»t)  “  or  -  F(“crt) 


Therefore,  for  K  >  0, 


F<5)  =  a 


For  t  >  0, 


(XjW.t)  -  op  -  0r  +  vp  -  F^OO  -  (vp+Cp)^ 


If  Vp  is  constant,  then  Xp(t)  =  vpt,  so  Cp  =  F^-Cptj.  Therefore,  for 
£  <  6,  F(5)  *  o^  +  Vp.  Summarizing, 


no 


a  if  5  >  0 
r 


Frcm  t  > 
that 


j  or  +  vp  if  ?  <  0 
0  and  from  the  formula  for  F(5)  where  ?  *  X  -  °-0r^,:>  ^°^OMS 


X  > 


(f^r) 


t4=>-x/t  >  c 


and 


X  < 


o-or^t4->-X/t  <  Cr  +  • 


Therefore, 


X  =  {^IT  °~°rj 

[X  *  i^T  a~° r) *J  A  [X  *  °'Ur)tj 


4=*^  2.X/t  >Cr  +  ^  vp 
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Hencer 

X/t  <  Cr  +  ^  vp  «=»o  =  or  +  vp 

and 

Cr  +  ^T1  VP  -  X/t  -  Cr  “►X/tfor  ^i  =  0 

and 

X/t  >  C  o*.a 
r  r 

Thus  a(X,t)  is  completelj'  defined.  Now  the  relation  between  o  and  p  is  one 
to  one.  so  p  is  completely  defined.  As  mentioned  previously,  P,  v,  c  are  all 
functions  of  p,  so  P,  v,  c  are  completely  defined.  Notice  that  °p  =  °r  +  vp» 
therefore,  op  will  ba  positive  if  jvpj  <  a Since  op  »  2/y-l  Cp,  for  the 
piston  not  to  leave  the  gas  (i.e.,  no  vacuum)  it  is  required  that  |vp|  <  2/y-l. 

In  the  left  side  constant  wave,  the  fluid  velocity  is  vp.  In  the  rarefaction 
wave,  the  fluid  velocity  goes  linearly  from  vp  to  Q.  In  the  right  side  constant 
wave,  the  fluid  velocity  is  zero.  See  figure  2c»  Since  there  is  a  one-to-one 
relation  between  velocity  and  density,  all  variables  may  be  written  in  terms 
of  the  velocity  because  by  the  equation  of  state  P  is  a  function  of  p, 

Xp(t)  »  XyCO)  +  vpt 

Xr(')  -  M0)  +  (Cr  +  ^ 

xc(t)  *  XpCO)  +  Crt 

For 

Xp(c)  <  X  <  XR(t) 
v(X,t)  «*  vp 
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X*(v)  = 


V3^  Vxc 


°-vp  vp 


So  the  kinetic  energy  fronX^  to  X^  is 


1/2 


+  HB3  +  5B,+|i] 


10B2 


where 


Y-l  VP 
B  «  — 

2  C 

r 


Then  the  total  internal  energy  is  the  total  energy  minus  the  total  kinetic 
energy. 

c.  Application  as  a  Test  ProblSc 
(1)  Eulerian  Input 

(a)  Initial  values 

For  0  <  X  £  XpCO) 


For  X  >  XpfO) 


v  =  0 


3c 
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P  =  P 


p  =  p 


(b)  Boundary  values 
For  X  =  0 


v  =  -  v. 


Y-l 


0  «  P^l  -  Ei  1^1  \ 


2 

y-i 


(2)  Lagrangian  input 
(a)  Initial  values 
For  x  >  0 


v  *  0 


P  -  P 


p  =  p 


(b)  Boundary  values 
For  x  *  0 


v(0,t)  =  -|vpj 

X(0,t)  =  X(0,0)  -|vp|t 
or,  since  X^t)  5  X<0,t),  then 


Xp(t)  =  XpCO)  -{vp|t 
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(3)  Numerical  Values  for  SCTP-II 

(a)  SCTP-II-A 

Pr  =  101*  dynes /cm2 

Pr  =  10”6  gm/cm3 

C2  =  y?  V  =  1.4  x  1010  cm2/sec2 
r  r  r 

|vp|  -  ',r/(7+ 1) 

AX  =  100  cm 

X_(0)  =  100  meters 
X^  =  300  meters 

A  reasonable  output  recipe  is  to  run  for  0.1  second  and  print  every  0.01  second. 
This  takes  about  120  cycles  for  a  CFL  number  of  1. 

(b)  SCTP-I.l-B 

Same  as  A  except  |vp|  ~  2Cr/(y+l) 

<c)  SCTP-II-C 

Same  as  A  except  Jvpj  -  2Cr/(y-l) 

(d)  SCTP-II* D 

Same  as  A  except  |vp|  =  4Cr/(y-l) 

(e)  SCTP-II-E 

Same  as  A  except  free  boundary  condition  on  the  left  in  place 
of  withdrawing  piston  on  the  left. 

(4)  Comments  on  the  Computer  Solution 
(a)  General 

If  specific  volume  is  used  in  a  code  instead  of  density,  then 
it  cannot  solve  those  problems  where  vacuums  occur  (C,  D,  and  E)  because  the 
specific  volume  becomes  infinite.  Whereas,  if  density  is  used,  the  vacuum 
'’ondltl  n  is  expressed  by  the  density  going  to  zero.  Also,  the  Eulerian 
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formulation  cannot  give  the  correct  answers  because  in  the  equation  of  motion 
the  density  (as  a  fun  on  of  time)  is  used  as  a  divisor.  (In  Lagrangian 
formulation  it  is  the  initial  density  which  is  used  as  a  divisor.)  The  mathe¬ 
matical  solution  to  C,  D,  and  E  is  the  same.  However,  it  cannot  be  expected 
that  the  computer  solutions  will  be  the  same.  Notice  that  in  C,  the  left 
boundary  of  the  gas  is  forced  to  move  at  the  correct  escape  velocity.  In  D, 
the  left  boundary  of  the  gas  is  constrained  to  move  too  fast  to  the  left.  In 
E,  the  code  is  allowed  to  calculate  what  the  escape  velocity  should  be.  Because 
the  masses  of  the  zones  of  the  fluid  do  not  go  to  zero  at  the  left  boundary, 
the  escape  velocity  that  the  code  computes  will  be  too  low.  That  is,  in  the 
mathematical  solution,  the  fluid  starts  moving  with  the  escape  velocity  instantly 
upon  its  release.  For  a  nonzero  mass  of  fluid  to  jump  from  zero  velocity  to  a 
constant  nonzero  velocity  requires  an  infinite  force.  The  mathematical  con¬ 
tinuum  model  is  not,  however,  predicting  such  a  thing.  The  continuum  model  is 
requiring  that  an  "infinitesimal  mass"  start  moving  with  a  constant  velocity 
and  is  thus  not  requiring  an  infin^-.e  force.  To  state  this  precisely,  let  x^ 
be  the  label  of  the  gas  point,  to  the  left  of  which  is  a  vacuum,  i.e., 

X^,t^  *  XR(t).  Let  xr  be  the  label  of  the  gas  point,  which  at  time  t  is  at 
the  boundary  of  the  undisturbed  gas  and  the  gas  which  has  started  moving  toward 
the  vacuum,  i.e.,  X , tj  »  X^,(t).  Let  x^  be  the  label  of  some  fixed  gas  point 
such  that  in  the  time  interval  of  interest  [0,T],  T  >  0,  the  gas  point  labeled 

x  is  quiet  or  undisturbed,  i.e. ,  x  >  x  for  0  <  t  <  T  (x„  is  a  function  of 
q  .  /  <3  r  -  -  r 

time),  e.g.,  let  X^x  ,fcj  *  X^,  Now  the  force  on  the  left  boundary  of  the  gas 
at  time  zero  is  the  negative  (the  force  is  directed  to  the  left)  of  the  initial 
pressure  P  of  the  gas  times  the  cross-sactional  area  of  the  pipe  (which  is  a 
unit  area)  and  this  is  equal  to  the  time  rate  of  change  of  the  momentum  of  the 
gas  at  time  zero.  Expressed  symbolically,  this  is 


-P  «*  lim 


XM 

d  f 


p(X,t)  v(X, t)dx 


-I  \ 

where  f^X(x,t),t)  5  f(x,t),  or  by  a  change  of  variables  in  the  integral,  we  have 
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-P  =  lim  jz  ^  p(k>0)  v(x, t)dx 

t-0  dt  J 
t>0  X£ 


(This  can  be  further  reduced  to 


r 

J  v(x •' 


-P  =  plim  - - - - 

t->0  c 

t>0 

by  defining  x^  *  0.)  So  that,  -although  the  time  derivative  of  v  at  x^  and  time 

zero  is  not  defined,  the  time  derivative  of  the  momentum  is  defined  and  finite. 

So  in  the  mathematical  continuum  model,  the  jumping  of  the  left  hand  boundary 
to  a  nonzero  constant  velocity  at  time  zero  does  not  require  an  infinite  force. 
However,  in  the  discrete  model,  the  finite  zoning  produces  a  positive  mass  for 
the  left  boundary  so  that  it  cannot  jump  from  a  zero  to  a  nonzero  constant 
velocity  by  the  application  of  a  finite  force.  The  finite  force  used  in  a 
code  should  be  just  the  negative  of  the  initial  pressure  at  the  left  hand 

boundary.  Then  the  acceleration  of  the  left  boundary  is  -P  divided  by  the  mass 

associated  by  the  code  with  the  left  zone  boundary  (which,  for  example,  in  PUFF, 
is  1/2  the  left  zone  mass).  In  C  and  D,  there  will  be  a  couple  of  compensating 
e  .’rcrs  influencing  the  total  energy  computation.  An  error  that  will  tend  to 
reduce  the  total  energy  is  caused  by  the  fact  that  the  pressure  at  the  left 
zone  boundary  Pp  will  not  be  zero;  therefore,  the  work  done  on  the  gas  at  the 
left  boundary,  Ppvpt,  will  be  negative.  This  causes  an  internal  energy  decrease 
by  the  first  law  of  thermodynamics.  An  error  that  will  tend  to  increase  the 
tota.  kinetic  energy  and,  therefore,  the  tocal  energy  is  the  numerical  integra¬ 
tion  of  the  kinetic  energy.  The  typical  coJe5s  numerical  integration  of 


1/2  j  p(X)  v2 (X)dX 


will  yield  too  large  a  value  for  the  kinetic  energy  because  this  will  be  approxi¬ 
mated  by  something  like 
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i*t*x 

1/2  v2(XX*1/2)  \ 


where 


XI+l/2  is 


XI+1  +  XI 


This  approximation  is  too  large  because  of  the  characteristics  of  the- particular 
p  and  v  encountered  in  this  problem.  This  error  affects  both  C  and  D.  A 
further  erroneous  increase  that  arises  in  D  is  caused  by  the  left  boundary's 
overly  high  velocity.  Thus,  the  kinetic:  energy  of  the  left  zone  is  erroneously 
increased;  therefore,  the  errors  affect  the  energy  partition  by  erroneously 
increasing  the  kinetic  and  decreasing  the  interal  energy.  The  total  energy 
may  either  increase  or  decrease,  depending  on  which  error  dominates.  Also, 
the  discrete  models  will  not  allow  the  density  to  go  to  zero  at  the  piston  face 
because  the  left  zone  started  with  a  finite  mass  and  is  just  being  stretched 
out  as  the  problem  progresses.  Other  points  of  interest  in  the  continuum 
solutions  which  should  be  approximated  by  the  discrete  solutions  are  these: 

In  A,  ^(t)  moves  to  the  right  with  velocity  Cf/2.  In  B,  XR(t)  is  stationary. 

In  C,  D,  E,  Xr(t)  moves  to  the  left  with  velocity  -2Cr/(y-l)  and  XR(t)  marks 
the  boundary  between  vacuisn  and  gas*  The  gas  can  just  keep  up  with  the  piston 
in  C,  i.e.,  Xp  -  X^,  and  at  the  piston  face,  the  values  of  the  pressure  and 
density  are  zero.  In  D,  the  piston  leaves  the  gas  behind  with  a  vacuum  between 
the  piston  at  Xp  and  the  gas  front  at  XR.  The  gas  front  can  move  into  a  vacuum 
with,  at  most,  the  velocity  -2Cr/(y-l)»  which  is  known  as  the  escape  velocity. 
So,  in  C,  D,  E,  Xp(t)  -  Xp(0)  -  2Crt/(Y~l)*  The  total  energy  from  the  piston 
to  an  undisturbed  point  on  the  right  increases  (or  decreases)  by  Ppvpt 
where  Pp  and  vp  are  the  pressure  and  velocity  at  the  piston  face.  In  A,  B, 

Pp  >  0,  but  in  C,  D,  E,  ?p  °  0,  so  in  A,  B,  the  total  energy  should  decrease, 
but  in  C,  D,  E,  it  should  remain  constant  from  start  to  finish. 

(b)  PUFF  comments 
1.  SCTP-1I-A 

The  typical  PUFF  profiles  in  velocity,  pressure,  and 
density  deviated  from  the  exact  solution  in  the  manner  described  in  figure  2d. 
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Pressure 

density  velocity 


Figure  2d.  PUFF  SCI.'  -II  Error  Graph 

Reading  figure  2d  from  right  to  left,  we  first  see  an  underround  at  Xc,  then 
an  over round  followed  by  an  undershoot  at  X^.  At  time  0.1  (170  PUFF  cycles) 
second,  the  velocity  deviations  at  the  undershoot  and  the  underround  are  both 
about  3  percent  (the  base  for  percentages  was  the  piston  velocity).  If  we 
take  the  initial  pressure  and  density  as  bases  for  the  percentage  deviation  in 
the  pressure  and  density,  respectively,  then  we  get  about  2  percent  as  the 
maximum  error  for  the  pressure  and  density.  The  total  energy  •.  imputations 
yielded  the  following: 

Xq  was  chosen  200  meters  to  the  right  of  Xp(0).  This 
yielded  a  total  internal  energy  at  time  zero  of  5  x  108  ergs.  Since  the 
velocity  and,  therefore,  the  total  energy  is  zero  at  time  zero,  the  initial 
total  energy  is  the  same  as  the  Initial  internal  energy.  After  170  cycles, 

PUFF  reached  the  problem  time  0.1  sec  and  the  energy  deviations  are  presented 
in  table  IJ-A.  Refer  to  the  energy  numbers  in  table  II-A.  These  numbers  are 
X108  ergs. 

PUFF  took  31  seconds  central  processor  time  on  the  CDC 
6600  and  170  cycles  to  run  this  problem  to  the  problem  time  of  0.1  second.  With 
a  CFL  number  of  1,  it  would  have  taken  120  cycles,  sc  that  PUFF's  effective 
CFL  is  120/170  ~0.7  on  this  problem. 
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2.  SCTP-II-B 

See  figure  2d  again  for  the  typical  PUFF  profiles.  The 
results  on  percentage  underround  and  undershoot  are  about  the  same  as  in  SCTP- 
II-A.  The  total  energy  computations  yield  the  following: 

Xq  was  again  chosen  200  meters  to  t  he  right  of  Xp(0); 
the  initial  total  internal  and  total  energy,  therefore,  are  5  x  108  erg3. 

After  170  cycles,  PUFF  reaches  the  problem  time  0.1  second  and  table  1I-B 
displays  the  energy  deviations. 

PUFF  took  19  seconds  CDC  6600  central  processor  time  on 
this  problem.  Again,  PUFF’s  effective  CFL  is  120/170  ~  0.7  on  this  problem. 

3.  SCTP-II-C 

See  table  II»C  for  the  error  numbers. 

This  took  PUFF  170  cycles  and  21  seconds  CDC  6600  computer 
time  to  tun  to  the  problem  time  of  0.1  second.  Thus,  the  effective  CFL  number 
is  120/170  ~  0.7, 

4.  SCTP-II-D 

See  table  II-D  for  the  error  numbers. 

This  problem  took  PUFF  170  cycles  and  30  seconds  CDC 
6600  computer  time  to  run  to  the  problem  time  of  0.1  second.  The  effective 
CFL  number  is  120/170  ~  0.7. 

5.  SCTP-II-E 

As  discussed  earlier,  the  gas  front  at  the  free  boundary 
of  the  gas  in  a  finite  difference  solution  will  not  move  instantly  with  the 
escape  velocity  because  the  left  zone  has  a  nonzero  mass.  This  causes  a 
"truncation  of  the  profile."  That  is,  the  graphs  appear  about  the  same,  but 
they  are  not  extended  as  far  to  the  left  as  the  exact  solution.  See  table  1I-E 
for  the  error  numbers. 

PUFF  took  170  cycles  and  29  seconds  CDC  computer  time  to 
run  to  0.1  second.  Thus,  effective  CFL  number  is  120/170  ~0.7. 
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Table  Il-A 


PUFF  ERROR  TABLE  FOR  SCTP-II-A 


Problem  time  *  0.1  sec  FUFF  cycle  =  170 

Computer  time  =  31  sec  Number  cf  active  zones  =  134 


Position  of 


Sum  Sqr.  Err. 

Max.  Err. 

Max.  Err. 

Pressure 

0,0032 

-0.0136 

xc 

Velocity 

0.0076 

+0.0248 

XR 

Density 

0.0025 

-0.0097 

xc 

Sum  Int.  Energy 

Sum  Kin.  Energy 

Sum  Tot.  Energy 

EXACT 

4.629  x  108 

1.027  x  107 

4.732  x  108 

PUFF 

4.629  x  108 

1.025  x  107 

4.731  x  108 

Table  II-8 

PUFF  ERROR  TABLE  FOR  SCTP-II-B 

Problem  time  «  0.1 

sec 

PUFF  cycle  -  170 

Computer  time  *  19 

sec 

Number  of  active  zones 

*  134 

Sum  Sqr.  Err. 

Max.  Err. 

Position  of 
Max.  Err. 

Pressure 

0.0032 

-0.0146 

xc 

Velocity 

0.0048 

-0.0178 

XR 

Density 

0.0027 

-0.0104 

xc 

Sum  Int,  Energy 

Sum  Kin.  Energy 

Sura  Tot.  Energy 

EXACT 

4.432  x  108 

2.923  x  107 

4.725  x  103 

PUFF 

4.431  x  10® 

2.920  x  107 

4.723  x  10° 
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Problem  time  ■  0.1  sec 
Computer  time  «  21  sec 


Table  II-C 

PUFF  ERROR  TABLE  FOR  SCTP-I1-C 

PUFF  cycle  <=  170 

Number  of  active  zones  =  135 


•Sum  Sqr,  Err. 

Max.  Err. 

Position  of 
Max.  Err. 

Pressure 

0.0036 

-0.0157 

XC 

Velocity 

0.0006 

-0.0022 

x 

Density 

0.0027 

-0.0112 

c 

xc 

Sum  Int.  Energy 

Sum  Kin.  Enerev 

Sum  Tot.  Energy 

EXACT 

4.260  x  108 

7.395  x  107 

5.000  x  108 

PUFF 

4.249  x  108 

8.045  x  107 

5.053  x  108 

Table  II-D 

PUFF  ERROR  TABLE  FOR  SCTP-II-D 
Problem  time  0.1  sec  PUFF  cycle  m  170 

Computer  time  -  30  sec  Number  of  active  2ones  „  U5 


Pressure 

Velocity 

Density 


EXACT 

PUFF 


Sum  Sqr.  Err. 

Max.  Err. 

Position  of 
Max.  Err. 

0.0038 

-0.0158 

X 

C 

0.0007 

-0.0023 

c 

0,0028 

-0.C113 

X. 

c 

Sum  Int.  Energy 

Sum  Kin.  Enertv 

4.260  x  108 

7.395  x  IQ7 

5.000  x  108 

4.245  x  10s 

9.919  x  107 

5.237  >:  108 
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Table  II-E 

PUFF  ERROR  TABLE  FOR  SCTP-II-E 
Problem  time  *  0.1  sec  PUFF  cvcle  =  170 


Computer  time  = 

29  sec 

Humber  of  active  zones 

=  134 

Sum  Sqr.  Err. 

Max.  Err. 

Position  of 
Max.  Err. 

Pressure 

0,0025 

-0.0140 

xc 

Velocity 

0.0006 

-0,0048 

free  boundary 

Density 

0.0019 

-0.0100 

xc 

Sum  Int.  Energy 

Sum  Kin.  Energy 

Sum  Tot.  Energy 

EXACT 

4.260  x  IQ8 

7.395  x  107 

5.000  x  108 

PUFF 

4.261  x  10a 

7.374  x  107 

4.998  x  108 
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3.  HYDROCODE  TEST  PROBLEM  SCTF-III 
a.  Description  of  Problem 

In  this  problem,  a  piston  proceeds  with  a  constant  acceleration  into  a 
gas  initially  at  rest  (by  "gas  initially  at  rest*'  is  meant  that  the  initial 
conditions  are:  velocity  is  zero;  density,  pressure  and  all  other  fluid 
parameters  are  constant).  This  forms  what  is  called  a  compression  wave.  At 
time  tg  =  2Cr/a(y+l),  a  shock  wave  is  formed  (tg  =  time  of  shock  formation, 

=  sound  speed  of  the  gas  at  rest,  a  =  acceleration  of  the  piston).  Until 
time  tg,  the  variables  are  continuous  and  the  solution  is  easily  found.*  More¬ 
over,  except  for  one  point  (the  front  of  the  compression  wave),  the  variables 
are  smooth  prior  to  t  .  Figure  3  attempts  a  graphical  representation  and 
description  of  the  problem.  As  shown  in  figure  3b,  the  compression  wave  front 
up  to  time  t  is  X£(t)  *  crt.  After  that  time,  the  compression  wave  front  is 
a  shock,  i„e.,  there  is  a  discontinuity  in  pressure,  density,  velocity,  etc. 


Compression  wave 


U.  v  *  at 
Xp  *  1/2  at2 

Piston  face 


Velocity  =  v^  *  0 

_  Dens it v  =  c  *  const. 

Gas  at  rest  Pressu;.e  «  pr  .  const. 

Sound  speed  =  =  const. 


X  =  C  t 

t 

Compression  wave  front 


Figure  3a.  Pipe  Plot  for  SCTP-III 


*See  K.  0.  Friedrichs’  paper  in  1948  C.P.A.M.,  page  211,  for  an  investigation 
of  the  solution  after  shock  formation. 
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Figure  3b,  Eulerian  Wave  Plot  for  SCTP-III 


b.  Derivation  of  Solution 


We  present  here  the  solution  for  the  preshock  region.  The  preshock 
region  is  the  set  of  all  points  (X,t)  such  that  0  _<  t  _<  tg,  i.e.,  the  set  of 
■points  in  the  X,t  plane  below  the  line  t  «  t  and  above  the  line  t  *  0.  This 
derivation  depends  on  three  key  points.  First  key:  there  is  a  one-to-one 
relation  between  the  density  and  the  velocity  and  the  pressure  of  the  fluid. 
Second  key:  the  values  of  the  density.,  pressure,  and  velocity  are  known  at 
the  piston  face.  Third  key:  the  density  surface  nay  be  described  in  terms 
of  its  level  lines  in  the  X,t  plane. 


Now  we  will  establish  the  first  key.  In  the  preshock  region,  v  -  o  >= 


v  -  o  *  -o  .  This  is  because  the  X  characteristics  (dr/dt  *  v  -  c)  extend 
r  r  r 

from  the  t  *  0  line,  where  v  *  v^  *  0  and  a  «  into  th  -•  preshock  region. 


and,  of  course,  v  -  o  is  constant  on  X_  characteristics.  Therefore,  the  pre¬ 
shock  region  is  a  X_  simple  wave,  i.e.,  v  *•  a  is  the  same  constant  in  the 
entire  region.  And  the  part  of  the  simple  wave  between  Xp(t)  and  X£(t)  is  a 
compression  wave  (Xp(t)  is  piston  face  position,  Xc(t)  is  sound  signal  front) 
Now,  for  0  <  t  <  t  ,  v(p)  ■=  o(p)  -  o  ,  o(p)  «  2c(p)/y-l,  c?(p)  «*  dP/dp  (p), 
c(p)  _>  0,  P(p)  °  kpY,  k  «  pr/PrY*  This  establishes  the  one-to-one  relations 
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mentioned  in  the  first  key  point.  Now  we  will  establish  the  second  key.  The 
velocity  at  the  piston  is  vp  -  at.  So,  o(p)  -  or  +  at  at  the  piston 
Then  utilizing  the  relations  between  o  and  c  and  ?  and  p,  the  second  h.-.y  is 
established.  So,  a-  the  point  ^Xp(t),t^,  i.e.,  the  piston  face,  we  know  the 
density,  pressure,  and  velocity.  By  the  first  key,  if  we  knew  p(X,t)  for  all 
values  ot  X  and  t,  we  would  know  the  complete  solution.  By  the  second  key,  we 
know  p (1/2  at2, t) .  The  third  key  suggests  that  we  represent  the  density 
surface  in  terms  of  its  level  lines.  Therefore,  we  investigate  the  level 
linas  of  the  density  surface.  That  is,  we  will  investigate  the  paths  in  the 
x,t  plane  on  which  the  density  is  constant.  To  be  precise,  we  should  denote 
points  on  this  path  fay  such  that  for  0  £  t  <_  tg  p^X^pc;t^,d  =  p 

That  is,  X  is  defined  to  be  a  function  of  two  variables,  density  and  time.  Or, 
as  some  phrase  it,  X  is  a  oue-parameter  (density)  family  of  functions  of  one 
variable  (time).  Let  us  denote  by  (x0,tQ)  one  point  such  that  p^tj  =  pq. 
So,  p(x(oo;t),tj  «=  p(x0,tQj  =  PQ  for  0  <  t  £  tg.  Now  the  level  line  Represen¬ 
tation  procedure  is  composed  of  two  parts.  First  part:  for  all  numbers,  p  , 
in  the  range  of  the  function  p,  find  numbers  X(p0)*t(po)  such  that  pfx/p 0),t(pQ' 
=  P0*  Second  part:  on  a  level  line,  ' 


so 


Therefore, 


dt 


.  dX 
°X  dt  + 


P=P0 


X  «  X 

o 


dt 


Then  let  Xq  =  X^p^,  c0  =  t(°c)v  Then  we  solve  for  pq  and  this  gives  us 
PQ(X,t)  which  Is  the  representation  of  the  density  surface  we  desire,  i-e., 
Po(X,t)  *  p(X,t).  At  this  point,  we  know  of  two  previous  points  which  might 
bother  the  reader  (we  hope  there  are  only  two).  First  point:  we  have  been 
saying  "level  line  representation"  and  to  be  precise,  perhaps  we  should  have 
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said  "level  curve  representation."  We  will  prove  that  the  level  curves  are 
lines.  Second  point: 


is  a  bit  vague  perhaps,  but  it  turns  out  that 


so  that 


d(vp) 

dp 


P=P 


o 


This  last  relation,  by  the  way,  proves  the  first  point  about  the  level  curves 
being  lines.  Now  we  will  establish  that 


d(vp) 

dp 


p*p 


PXP, 


In  a  smooth  flow,  conservation  of  mass  is  expressed  by  p  +  (pv)v  =  0.  By  the 

C  A 

first  key,  v  =  v(p),  so 


P 


t 


d(pv)  , 
dp  "X 


=  0 


Therefore, 

Pt  „  d (pv) 
'  px  ~dp 
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Therefore, 


_  d(py) 
dp 


~  d (pv) 
dt  “  dp  ' 


13  the  dIf£erenUal  e^ti0"  '«  “•  l-el  „uh  density  „ence_ 


‘M  ■ 


d(vp) 

dp 


,.p  ]  K) 

0_ 


dv  _  dv 

dp  V  +  p  do 


pH^=Cbyv  =  o- 


r>2  _  dP 
C  ~d£ 


P  =  kpY 


—  =  d(vo) 
dc  “  dp 


=  (v+c; 


/ 


i 
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Therefore,  the  level  li.ies  are  characteristics.  Thus,  in  short  summary, 
the  density  lavel  curves  are  lines  and  also  they  ara  A+  characteristics 
(figure  3c). 


A_  characteristics  have  slope  -C^. 

A+  characteristics  have  slope  (y+l)/2  at  +  at  the  point  (1/2  at2,t). 

Figure  3c.  Characteristic  Linas  Plot 

If  the  mininur*  intersection  time  of  the  A+  characteristics  is  computed,  it  is 
found  to  be  2Cr/a(y+l).  A  point  where  A+  characteristics  cross  may  indicate  a 
discontinuity.  This  is  because  v  +  a  is  a  different  constant  on  different  A+ 
characteristics.  So,  when  two  characteristics  cross,  the  contradiction  can 
sometimes  be  resolved  by  a  discontinuity  at  the  point.  And,  of  course,  the 
discontinuity  is  known  as  a  shock.  Now  we  have  X (p ; t )  and  the  problem  is  to 
find  p(X, t);  so  we  must  "invert"  the  relation 

x  ■  x° +  (v(°°) +  c(p°))  K) 
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v;here 


(X«'f'o) 


%  =  p 

By  the  second  key,  v,  c,  c,  p,  p  are  known  at  the  piston  face.  All  values  of 

the  velocity  that  are  taken  on  will  be  taken  on  at  some  time  at  the  piston  face, 

Therefore,  by  the  first  key,  all  values  of  o,  C,  P,  p  will  be  taken  it  some 

time  at  the  piston  face.  By  the  first  key,  if  we  know  the  velocity  surface, 

we  know  the  complete  solution.  It  is  computationally  convenient  to  switch  to 

the  velocity  as  the  function  to  determine  first.  By  the  first  key,  the  density 

level  lines  are  also  velocity  level  lines.  Thus,  for  any  value  of  the  velocity 

that  is  taken  on  in  the  preshock  region,  we  can  find  an  X  ,t  on  the  piston 

o  o 

x0»t0 J  -  vq.  Thus,  we  can  produce  the  x/v  t|vo}  such  that 

„l  V  \  .1,.  \  1  .  ..  ..LJ.L  i.  _ f _ ,  ,  .  .  O  l  .  '  1 


'(X(Vo)-C(Vo))  ■ 

representation  procedure. 


vq  which  is  required  by  the  first  part  of  the  level  curve 


X  *  1/2  at2  + 
o 


(“. +  c  (1/2  <•<))  ho) 


and 


t.  =  t(v  )  =  — 
o  \  oj  n 


Therefore, 


*1/2a(r)  +  (’.  +  c  (’'.))  («-r) 


Solving  for  v  ,  we  get 


vq  =  v(X, t)  = 


~(cr  -  1/2  (t+1)  at)  +  J  (cr  -  at)2  +  2ay  (crt-x) 


for 


and 


X  <  X  <  C  t 
P  --  -  r 


0  <  t  <  t 
—  —  s 
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Also,  v(X,t)  =  0  for  X  >  C  t  and  0  <  t  <  t  .  Then  using  the  simple  wave 

r  s 

formulas: 


C  =  C 


(**n) 


P  *  P. 


ft) 


P  «  p 


(f) 

\Cr) 


Notice  that  at  time  t  *  2C  /a(y+l)  and  position  X  =  C  t  ,  the 

8  IT  8  V  8 


lim  v, 


:K)  * 


This  indicates  that  a  shock  forms  at  ^Xg,tg^.  Now  we  compute  the  total  energy 
and  energy  partition.  The  total  energy  should  increase  according  to  the  amount 
of  work  done  on  the  gas  by  the  piston,  which  is 


Pp ( T ) Vp ( T ) d  T 


Integrating  by  parts,  this  becomes 


2c  P 

+ _ _ _ 

+  a(3y~l)(2y-l) 


AFWL-TP.-67-127 


To  get  the  energy  partition,  we  ha\>e  the  choice  of  computing  either  the  kinetic 
or  the  internal  and  subtracting  from  the  total  energy  to  gtt  the  other.  We 
will  compute  the  internal  energy  in  «:he  compression  wave. 


Internal  energy 


P 

Y“1 


dX 


by  the  equation  of  state.  Using  the  velocity  level  Line  equation  which  is  a 
quadratic  in  the  velocity,  we  get 


dX  = 


=  "  a  (cr  ~  at  +  YVj  dv 
At  X^,  v  ■=  v^  at.  At  Xc,  v  =  vr  =  0.  Using  these  substitutions,  we  get 


where 


Internal  energy 


P(v)(A+yv)  dv 


A  = 


c 


r 


at 


Then,  using  the  simple  wave  formula  for  P(v),  the  integral  is 


1 

a(y-i) 


(A+yv)  dv 


and  this  can  be  integrated  by  parts. 
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Internal  energy  = 


a(y-l) 


TvT 


&) 


(A+yv) 


(1  +  1rr) 

Mf^2] 


+  2 

Y-l 


(&)■ 


Finally,  kinetic  energy  *  total  energy  -  internal  energy.  For  computation 
purposes,  it  seems  best  to  choose  a  point,  X^,  far  enough  to  the  right  of  the 
piston  so  that  the  fluid  is  quiet  (at  rest)  at  X^  throughout  the  course  of  the 
problem.  Then  take  the  initial  total  energy,  which  will  be 


(W0))  pr  er  “  r^i  (W<0>) 


and  add  the  total  energy  increase, 


PDv  dx 
P  p 


This  gives  the  total  energy  from  Xp  to  X^.  Then  the  internal  energy  will  be 
the  Internal  energy  in  the  compression  wave,  computed  above,  plus  the  internal 
energy  from  X„  to  Xn,  which  is 

V/  V 


(VXc)  •  »r  •  er  ■  ( VXc)  7^ 
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(3)  Numerical  Values  for  SCTP-1II 

(a)  SCTP-I1I-A 

Pf  =*  lO4  dynes/cm2 

Pr  =  10-6  gm/cm3 

c2  ■=  yP  V  =  1.4  x  1010  cm2/sec2 
r  r  r 

a  «  c  /I  sec 
r 

Ax  =  10  meters 

Right  boundary  at  1500  meters 

A  reasonable  recipe  for  output  is  as  follows:  Run  for  1  sec  (problem  time  not 
computer  time)  and  print  at  0.1  sec  intervals  and  at  time  2/(y+l)  (when  the 
shock  forms).  If  the  CFL  number  equals  1,  this  takes  218  cycles.  This  follows 
from 


At 


n 


*  min 
I 


AX“  „ 
piston 

n 

cpiston 


Then  by  conservation  of  mass, 


AX 

Ax 


P 


So, 


At 


n 


By  the  simple  wave  formulas 


2 


and 
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Density 

Pressure 

Velocity 


Figure  3d,  PUFF  SCTP-III  Error  Graph 

1.  /-.SCTP-III-A 

ACCURACY:  See  table  III-A. 

TIMING:  PUFF  took  245  cycles  and  19  seconds  CP  time 
(CDC  6600)  to  run  to  0.8333  second.  Effective  CFL  number  =  218/324  «  0.68 
on  run  to  1.0  second. 

2.  SCTP-III-B 

ACCURACY:  See  table  III-B. 

TIMING:  PUFF  took  245  cycles  and  21  seconds  CP  time 
(CDC  6600)  to  run  to  0.08333  second.  Effective  CFL  number  =  218/324  «  0.68 
on  run  to  0.1  second. 
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Table  7.IX-A 

PUFF  ERROR  TABLE  FOR  SCTP-III-A 

Problem  time  = 

0.8333  sec 

PUFF  cycle  =  245 

Computer  time 

-  19  sec 

Number  of  active  zones 

=  19 

Sum  Scjr.  Err. 

Max.  Err. 

Position  of 

May 

Pressure 

0.0057 

0.0426 

iiOA  •  UL  L  • 

X^ 

Velocity- 

0.0125 

0.0909 

c 

x 

Density 

0.0053 

0.0407 

c 

xc 

Sum  Int.  Energy 

Sum  Kin,  Enerpv 

Sum  Tot.  Energy 

EXACT 

4.36800  x  109 

2.62864  x  108 

4.63087  x  109 

PUFF 

4.37190  x  109 

2.64055  x  108 

4.63595  x  109 

Table  1II-B 

PUFF  ERROR  TABLE  FOR  SCTP-III-B 


Problem  time  ~ 

0.08333  sec 

PUFF  cycle  =  245 

Computer  time  = 

:  21  sec 

Number  of  active  zones 

-  119 

Sum  Sqr.  Err. 

Max.  Err. 

Position  of 
MflX •  Err  - 

Pressure 

0.0057 

0,0426 

X_ 

Velocity 

0.0125 

0.0909 

C 

X. 

Density 

0.0053 

0.0407 

C 

xc 

Sum  Int.  Energy 

Sum  Kin.  Energy 

Sum  Tot.  Energy 

EXACT 

4,36800  x  108 

2.62864  x  IQ7 

4.63087  x  108 

PUFF 

4.37190  x  10° 

2.64055  x  107 

4.63595  x  108 

l 

* 

i 


61 


Figure  4b.  Wave  Plot  for  SCTP-IV 
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b.  Derivation  of  Solution 


The  analysis  is  the  same  as  for  SCTP-III. 


v(X,t) 


-(cr  -  at)  + 


!cr  -  at 


t- 

)  +  2aY( 


for  at2/2  <_  X  <_  c^t.  v(X,  t)  =  0  for  X  >_  cr.t.  For 


2c 


C  as  c  + 

r  2  v 


P  =  P 


« 


Y-l 


P  *  P 


4) 


2 

Y-l 


For 


v  > 


2c 
_ r 

Y-l 


there  is  a  vacuum;  therefore,  0  =  P  =  p  =  c  from  Xp(t)  to 


xp(tv)  -  (t-tv) 


t  is  the  time  the  vacuum  is  formed, 
v 


2c 

Sa|t  =  — f 
1  1  v  y-l 
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(2)  Lag rang ian  Input 

(a)  Initial  values 


P  =  P,p=p.v  =  v 
r»  *  r 

(b)  Boundary  values 
At  x  =  0, 


0 


X  =  -  1/2 ] a | t2 ,  v  =  -ja|t 

(3)  Numerical  Values  for  SCTP-IV 
(a)  SCTP-IV-A 

Pr  *  104  dynes/cm2 

p^  =  10" 6  gm/cm3 

c2  =  yp  V  »  1,4  x  1010  cm2/sec2 
r  '  r  r 

a  =  -  c  /I  sec 
r 

Ax  =  1  meter 

Right  boundary  at  150  meters. 

A  reasonable  output  recipe  is  to  run  for  10  seconds  with  prints  at  1,  4,  5,  6, 
and  10  seconds. 


(b)  SCTP-IV-B 

Same  as  A  but  a  =  -  10  cr/l  second  and  print  at  0.1,  0.4, 
0.5,  0.6,  and  1.0  second. 

(4)  Comments  on  the  Computer  Solution 

(a)  General  comments 

After  time 
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v  a 


there  is  a  vacuum  between  the  piston  and  the  gas  at 


2c  . 

1/2  at2  +  — £  (t-t  ) 
v  y-1  V  v/ 


the  pressure  and  density  are  zero  for  t  ty.  From  the  left  edge  of  the  gas  at 


1/2  at2  +  — f  (t-t  ) 

v  y-1  \  v/ 


the  same  formulas  hold  as  before  tv>  If  the  hydrocode  uses  specific  volume  as 
a  variable  instead  of  density,  there  will  be  an  overflow  when  the  density  goes 
to  zero.  The  total  energy  should  decrease  until  the  vacuum  is  formed,  but 
after  t^,  the  total  energy  should  not  change.  The  energy  decrease  Is  the  work 
done  at  the  piston 


U 

f  Pp(t)vp(t)dt 


This  was  computed  in  the  derivation  of  solution.  A  Lagrangian  code  with  zone 
center  pressures  will  not  compute  the  energy  decrease  correctly  because  the 
piston  face  zone  will  never  have  zero  pressure.  The  energy  partition  should 
also  be  checked  against  the  EXACT  solution  as  computed  in  the  derivation  of 
solution.  For  further  comments,  look  back  to  the  comments  in  SCTP-II. 

(b)  PUFF  comments 

As  previously  observed,  it  is  impossible  for  PUFF  to  conserve 
energy  on  this  problem.  The  worst  behaving  variable  is  the  nonprimary*  variable, 
momentum;  the  reason  this  is  off  so  much  is  because  its  computation  involves 


*Primary  and  nonprimary  variables  are  defined  In  Appendix  III,  Error  Functions 
and  Error  Formats. 
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the  mixing  of  a  zone  center  quantity  (density)  and  a  zone  boundary  quantity 
(velocity),  thus  introducing  interpolation  errors  which  cannot  be  ascribed  to 
PUFF.  PUFF's  primary*  variables  are  quite  near  the  EXACT  solution.  See 
figure  4c  and  tables  IV-A  and  IV-B  for  accuracy. 

TIMING:  PUFF  took  25  seconds  CDC  6600  CP  time  and  169  cycles 
to  run  to  10.0  seconds  on  SCTP-IV-A.  On  SCTP-IV-B,  PUFF  took  16  seconds  CDC 
6600  CP  time  and  19  cycles  to  run  to  1.0  second. 


PUFF  Primary  Variables: 


Figure  4c.  PUFF  SCTP-IV  Error  Graph 


*See  Appendix  III  for  definition  of  primary  variables. 
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Table  IV-A 


PUFF  ERROR  TABLE  FOR  SCTP-IV-A 


Problem  time  =  10  sec  PUFF  cycle  =  169 

Computer  time  =  25  sec  Number  of  active  zones  =  134 


Pressure 

Velocity 

Density 


Position  of 


Sum  Sqr.  Err. 

Max.  Err. 

Max.  ! 

0.0017 

-0.0099 

xc 

0.0004 

-0.0038 

xr 

0.0013 

-0,0071 

EXACT 

PUFF 


Sum  Int.  Energy 
3.09266  x  1010 
3.04551  x  1010 


Sum  Kin.  Energy 
6.98426  x  109 
8.92017  x  109 


Sum  Tot.  Energy 
3.79108  x  1010 
3.93752  x  1010 


Table  J.V-B 


PUFF  ERROR  TABLE  FOR  SCTP-IV-B 

Problem  time  *  1  sec  PUFF  cycle  =  19 

Computer  time  *  16  sec  Number  of  active  zones  =  21 

Position  of 


Sum  Sqr.  Err 

Max.  Err. 

Max,  Err. 

Pressure 

0.0166 

-0.0492 

xc 

Velocity 

0.0018 

-0.0042 

xc 

Density 

0.0132 

-0.0358 

xc 

Sum  Int.  Energy 

Sum  Kin.  Energy 

Sum  Tot .  Energv 

EXACT 

3.684  x  10l° 

6.984  x  103 

3.754  x  1C10 

PUFF 

3.671  x  1010 

2.844  x  10s 

3.956  x  10i0 
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5.  hfDROCODE  TEST  PROBLEM  SCTP-V 


a.  Description  of  Problem 


This  is  called  the  shock  tube  problem.  It  an  example  of  the  more 

general  Riemann  problem.  The  Riemann  problem  is  that  of  determining  the  flow 
after  the  conjunction  of  two  states,  left  state  and  right  state,  with  ?£,  p„, 
the  constant  values  of  the  left  state  and  P  ,  p.r,  the  constant  values  of 
the  right  state.  In  tne  shock  tube  problem,  and  v^  are  no  longer  arbitrary 
but  are  set  to  zero.  So  the  problem  may  be  interpreted  as  the  determination 
of  the  flow  after  removal  of  a  membrane  separating  two  constant  states  at  rest. 
We  choose  to  make  as  a  convention  P.  >  P  >0.  Then  in  the  code  test  problem 
we  will  run  through  the  three  possibilities 


p£  >  pr 


p£  <  °r 


At  time  zero,  the  membrane  is  removed.  The  resultant  action  is  a  rarefaction 
wave  traveling  into  the  left  state  and  a  shock  traveling  into  the  right  state. 
The  velocity  is  v£  *  G  to  the  left  of  the  rarefaction  wave.  From  the  left  of 

the  rarefaction  wave  to  the  right,  the  velocity  rises  linearly  from  0  to  v  >  0. 

The  velocity  is  constant  at  v  from  the  right  of  the  rarefaction  wave  rightward 
toward  the  shock.  At  the  shock,  the  velocity  jumps  from  v  >  0  down  to  v^  =  G. 
See  figure  5d  for  a  velocity  plot.  The  pressure  drops  continuously  across  the 
rarefaction  wave  from  P.  to  P  .  The  pressure  has  the  value  P  constantly  from 
the  right  of  the  rarefaction  wave  to  the  shock.  The  pressure  drops  from  P^  to 
P^  across  the  shock*  See  figure  5e  foi  a  pressure  plot.  The  density  drops 
continuously  across  the  rarefaction  wave  from  p,  to  a  o  ,  which  value  it  main- 
tains  from  the  right  of  the  rarefaction  wave  to  the  point  in  the  fluid  where 
the  initial  discontinuity*  was  and  there  the  density  jumps  up  to  c  .  From  the 

initial  discontinuity  point  to  the  shock,  the  density  is  p  .  At  the  shock. 


*This  density  discontinuity  is  called  a  contact  discontinuity. 
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Membrane 

l _ 


Figure  5a.  Pipe  Plot  of  Initial  Conditions  for  SCTP— V 


Figure  5b.  Eulerian  Wave  Plot  for  SCTP-V 
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t 

„Ra re faction 

H  wave  j  ^ 

Contact 

^.discontinuity 

path 

\  \ 

vP-Shock  path 

Lagrangian  Wave  Hot 


Figure  5d.  Velocity  Plot 
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Figure  5f.  Density  Plot 
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the  density  jumps  down  from  to  p^,  which  value  is  maintained  all  the  way  to 
the  right.  See  figure  5f  for  a  density  plot.  We  will  mention  again  that  this 
is  a  special  case  of  the  Riemann  problem  because  in  the  Riemann  problem,  the 
velocities  are  also  arbitrary  along  with  Pc«  p0,  P  ,  p  .  Figures  5a  through 
5i  give  graphic  explanations. 

b.  Derivation  of  Solution 

Assume  states  m  and  r  are  connected  by  a  shock  facing  to  the  right,  i.e.. 
?m  >  F  .  By  the  conservation  of  momentum  relation,  (2:N2L),  derived  for  SCTP-I, 


=  P  -  P 
m  r 


m  (v  -v  i  -• 

V  m  x) 

By  the  conservation  of  mass  relation,  (1:CM),  derived  for  SCTP-I, 


v 

r 


v 

s 


m 


v 


m 


v 


s 


m 


So 


00 


w2(V  -V  )  =  P  -  ? 

\  r  m  /  m  r 
and 


m  V  -V 
r  m 

By  the  Rankine-Hugoniot  relation,  (5:RH),  derived  for  SCTP-I,  the  may  be 
eliminated  to  yield 


! 

I 


m 


F  +  P  (y-M)  P  +  (y-1)  P 

2  =  m  r  y+1  _  m  1  _ r 


y+1  ’7r 


ZV 
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\ 
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So,  for  a  rarefaction  wave  moving  to  the  left. 


We  summarize  and  review  the  action:  In  the  shock  tube  problem,  if  Pfc  >  Pr>  a 
shock  will  travel  to  the  right  and  a  rarefaction  wave  will  travel  to  the  left. 
To  the  right  of  the  rarefaction  and  the  left  of  the  shock  will  be  a  region 
constant  in  pressure  and  velocity  but  having  one  discontinuity  in  density  at 
the  point  in  the  fluid  of  the  initial  discontinuity.  Other  than  this  jump, 
the  density  is  constant  (it  takes  on  only  two  values)  in  the  middle  region. 

See  figures  5d,  5e,  and  5f.  Proceeding  from  left  to  right,  the  density  is  p£ 
in  the  left  constant  wave;  then  it  drops  continuously  throught  the  rarefaction 
wave  to  p  then  it  jumps  up  across  the  initial  discontinuity  to  p  J  then 
it  jumps  down  as  it  crosses  the  shock  to  p^.  If  the  initial  discontinuity  has 
the  Lagrangian  label  x  -  0,  then  the  line  x  =  0  in  the  Lagrangian  wave  plot 
is  the  path  of  the  initial  discontinuity.  Therefore,  the  line  x  =  0  is  the 
path  of  the  contact  or  tangential  discontinuity.  See  figure  5c  for  the 
Lagrangian  wave  plot. 

Shock  relation: 

v  *  v  +  ♦  (p  ) 
m  r  r\  m) 
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Rarefaction  relation: 


v  ”  vt  -  V0(?  ) 
m  l  £.  \  m/ 

These  two  relations  determine  v  and  P  .  The  equation  for  P  obtained  by 

mm  m 

equating  the  two  right  sides  of  the  above  relations  is 


2V 


1 

2y 


(y+D  Pm  +  (y-1)  Pr 


*  + 


/»» 

z 


-  p 


m 


or 


or  for  y  =  1*4> 


ml 


is  determined  by 


the  isentropy  condition 


P 

m 


V  =  p  1  is  determined  by  (5:RH): 
mr  mr 


0  =  e 
mr 


77 


AFWL-TR-67-127 


where 


The  shock  velocity  is  determined  by  the  shock  speed-pressure  relation,  (7:SS(P))> 
derived  for  SCTP-I: 


The  velocity  in  the  rarefaction  wave  linearly  connects  v  =  0  to  v  as  was 

&  ISl 

shown  in  SCTP-II.  Then  in  the  rarefaction  wave. 


2y 


as  was  shown  for  SCTP-II.  The  left  side  of  the  rarefaction  wave  is  at 

xc(t)  =  Xs(0)  -  C^t 

and  the  right  side  is  at 

=  V°>  -  (c4  -  *r 

as  was  shown  for  SCTP-II.  The  shock  wave  is  at 

X$(t)  =  Xg(0)  +  Vgt 

The  initial  discontinuity  point  of  the  fluid  is  at 

V0  3  Xs(0)  +  vat 
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Solution  summary: 


Region 


For  X  <_  Xc(t),  the  values  are  P  ,  p  v0 
^  x  x  x 


Xg(c)  £  ^  i  X^(t),  the  values  are 
X  -  X  (t) 

V(X,t)  =  X^(tj  -  Xc(t)  Vm 


Rarefaction 

Region 


c-c‘(1'^) 


P  -  P, 


ft)' 


Middle 

Region 


p  =  p 


ft) 


P°r  X_  <  X  <  Xc(t),  the  values  are  P  ,  v 
«•  o  m  tn 

For  X^(t)  £X  <  XjjCt),  the  density  is  p^ 
F°r  Xjj(t)  <  X  <  Xg(t),  the  density  is  p 


Right 

Region 


For  X  >  Xc(t),  the  values  are  P  ,  p  ,  v 
S  r*  r*  r 


c*  Application  as  a  Test  Problem 
(1)  Eulerian  Input 

(a)  Initial  values 


For  X  <  Xs(0), 
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fv 


i** 

i<; 


(b)  Boundary  values 

At  x  =  0,  hold  the  values  at  P.,  p„,  v. 

(3)  Numerical  Values  for  SCTP-V 

(a)  SCTP-V-A 

Xg(0)  =  100  meters 
AX  *  1  meter 
P^  =  10?  dynes/cm2 
P£  =  10“ 5  gm/cm3 

Vi  =  0 

P^  =  101*  dynes /cm2 

Pr  =  10“ 6  gm/cm3 

v  =  0 
r 

Right  boundary  at  250  meters 

These  imply  the  following  values: 

P  =  1.888  x  107  dynes/cm2 
o 

v  *  3.964  x  106  cm/sec 
m 

=  3.040  x  10“6  gm/cm J 

p  =  5.S82  x  1 0-6  gm/cm3 
mr 

.  4.760  x  106  cm/sec 

v  = 

s 

The  output  recipe  we  used  is:  print  out  at  time  1  x  10-*4  sec,  1  x  10' 
and  2  x  10“3  sec. 

(b)  SCTP-V-B 

Sane  as  A  except 
X_(0)  =  250  meters 
p^  =  10" 6  gm/cn3 
Right  boundary  at  500  meters 


These  imply  the  following  values: 
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=  4.610  x  10 7  dynes /cm2 

v  =  6.196  x  1G1  cm/sec 
□ 

=  5.751  x  10“  7  gm/cm3 
Pfflr  =  5.992  x  10~6  gm/cm3 
vg  ^  7.437  x  10^  cm/see 

(c)  SCTP-V-C 

Same  as  A  except 
Xs(0)  =  250  meters 
P£  =  10"6  gm/cm3 
Pr  =  10~5  go/cm3 
Right  boundary  at  500  meters 
These  imply  the  following  values: 

P^  =  7.406  x  ID7  dynes/cm2 
vq  *  2.484  x  XGS  cm/sec 
Pmi  8.070  x  10"  7  gm/cm3 

pmr  *  5*995  x  10~5  ga/c®3 
Vg  =  2.931  x  106  co/sec 

(4)  Comments  on  the  Computer  Solution 
(a)  General 


The  variations  A,  B,  C  were  made  to  explore  the  three 

possibilities  p,  >p,p,  =  p,p_<p.  The  solution  is  good  until  the 
v  r  i  r  t  r 

rarefaction  wave  reaches  the  left  hand  boundary  or  the  right  hand  boundary. 

At  the  last  time,  we  have  the  active  zone  situation  as  presented  in  table  V-A. 

(b)  PUFF 

ACCURACY:  On  SCTP-V-A,  the  most  noticeable  error  was  a 
smearing  of  the  density  discontinuity  at  Z^  (see  figure  5h).  The  only  other 
errors  were  the  typical  round  uaders  and  round  overs  at  corners.  (See 
figure  5h  and  table  V-B.  On  SCTP-V-B,  there  was  a  bit  of  oscillation  in  the 


av- 
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TIMING:  SCTP-V-A  took  976  cycles  and  61  seconds  CDC  6600 
CP  time  to  run  to  2  x  10" 3  sec.  SCTP-P-B  took  1527  cycles  and  168  seconds  CDC 
6600  CP  time  to  run  to  2  x  10“3  sec.  SCTP-V-C  took  611  cycles  and  71  seconds 
CDC  660(5  CP  time  to  run  to  2  x  10“ 3  sec.. 
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Table  V-A 

ACTIVITY  TABLE  AT  TIME  2  x  10“3  SECOND  ON  SCTP-V 


(meters)  (meters) 

SCTP-V-A  25  195 

SCTP-V-B  14  400 

SCTP-V-C  3.4  310 


Table  V-B 

PUFF  ERROR  TABLE  FOR  SCTP-V-A 


Problem  time  =  2  x 

10“ 3  sac 

PUFF  cycle  =576 

Computer  time  =  61 

sec 

Number  of  active,  zones 

=  202 

Sum  Sqr.  Err. 

Max.  Err. 

Position  of 
Max.  Err. 

Pressure 

0.0061 

-0.0545 

X 

s 

Velocity 

0.0795 

+0.875 

X 

s 

Density 

0.0271 

+0.172 

X 

s 

Sum  Int.  Energy 

Sum  Kin.  Energy 

Sum  Tot.  Energy 

EXACT 

2.177  x  1012 

3.222  x  10n 

2.500  x  1012 

PUFF 

2.178  x  1012 

3.218  x  101 1 

2.500  x  1012 
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Table  V-C 

TUFF  ERROR  TABLE  FOR  SCTP-V-B 

Problem  Lime  =  2  x  10“3  sec  PUFF  cycle  =  1527 

Computer  time  *.  168  sec  Number  of  active  zones  «  405 


Sum  Sqr.  Err. 

Max.  Err. 

Position  of 
Max.  Err. 

Pressure 

0.0137 

-0.263 

X 

s 

Velocity 

0.4447 

+0.748 

X 

3 

Densit}* 

0.0306 

-0.415 

X 

s 

Sum  int .  Energy 

Sum  Kin.  Energy 

Sinn  Tot.  Energy 

EXACT 

5.668  x  1012 

5.834  x  10Ji 

6.251  x  J.012 

PUFF 

5.669  x  1012 

5.814  x  10n 

6.250  x  1012 

Table  V-D 

PUFF  ERROR  TABLE  FOR  SCTP-V-C 

Problem  time  *  2  x 

10'3  see 

PUFF  cycle  -  611 

Computer  time  »  71 

sec 

Number  of  active  zones 

*  316 

Sum  Ser,  Err. 

Max.  Err. 

Position  of 
Max.  Err. 

Pressure 

0.0282 

-0.478 

X 

s 

Velocity 

0.0487 

+0.701 

X 

s 

Density 

0.0360 

-0.471 

X 

8 

Sum  Int.  Energy 

Sum  Kin.  Energv 

Sum  Tot .  Energv 

EXACT 

6.005  x  1012 

2.460  x  1011 

6.251  x  1012 

PUFF 

6.007  x  10i2 

2.434  x  1011 

6,250  x  1012 
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6.  HYDROCODE  TEST  PROBLEM  SCTP-VI 


a,  Description  of  Problem 


This  problem  is  the  collision  of  two  shock  waves.  Proceeding  from  left 

to  right,  the  values  are  P0,  p0,  v.  connected  by  a  right  facing  shock  to  P  , 

Pq,  v  which  in  turn  is  connected  by  a  left  facing  shock  to  P^,  p^,  vf.  As  a 

convenient  convention,  always  take  P  _>  P  .  Figure  6  gives  a  graphical 

36  IT 

description. 

b.  Derivation  of  Solution 

After  collision,  the  pressure  profile  looks  like  figure  6d.  Now  by 
the  analysis  in  SCTP-V,  v^  =  v^  +  for  a  shock  facing  to  the  right  and 

similarly,  vm  15  va  “  f°r  a  shock  facing  to  the  left,  where 


. 

♦a<»  3  (P-Pa)  ft+i)  P  +Vl)  P- 


From  these  two  equations,  P  and  v  are  determined.  Then  by  the  Rankine- 
’  m  m 

Hugoniot  relation,  (5:RH),  pm£  and  Pmr  are  determined.  Let 


v  =  p~J 


V  =  p  1 
mr  mr 


Then  by  (5:RH), 


Fl  (ftVVl m)  '  -V1 


~  (p  V  -F  V  \  =  — (v  -v  ) 

(- 1  \  m  mr  m  mr/  2  \  r  mr/ 


Let  t  be  the  time  of  collision,  X  (t)  be  the  position  of  the  left  shock 

C  S  36 

prior  to  collision,  i.e.,  t  <  t£,  then 
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Figure  6a.  Pipe  Plot  for  SCTP-VI  Initial  Values 


Figure  6b.  Lagrangian  Wave  Plot  for  SCTP-VI 


Figure  6c.  Pressure-Density  Plot  Prior  to  Collision 
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m 


Figure  6d.  Pressure  Plot  after  Collision 


Xs*(t)  *  XsH<0)  +  lV 


Likewise,  let  X  (t)  be  the  position  of  the  right  shock  for  t  <  t  ,  then 

SC  0 

x  (t)  =  x  (0)  +  v  t 

sr  sr  sr 

Let  X*^(t)  be  the  position  of  the  left  hand  shock  after  collision  ^t  >  t^J, 
then 


X*  (t)  =  X  ,  +  v*0 (t-t  ) 
s£  col  s£\  cl 


where  Xcq^  is  the  place  the  collision  occurs.  Likewise,  let  X*f(t)  be  the 
position  of  the  right  hand  shock  after  collision;  then 


X*  (t)  =  X  ,  +  v*  (t-t  ) 
sr  col  sr\  c/ 


For  t  <  t  , 
c 


X  <  Xg  (t)  implies  the  values  are  P^,  p^, 

X  .  <  X  <  X  (t)  implies  the  values  are  P  ,  p  ,  v 
s£  srv  Y  o’  *o*  o 

X  (t)  <  X  implies  the  values  are  P  ,  p  ,  v 
sr  v  r’  r’  r 


For  t  >  C  , 
c 


X* 


s£ 


X  <  X*  (t)  implies  the  values  are  P^,  p£,  v^ 

(t)  <  X  <  X  -  +  v  (t-t  J  implies  the  values  are  P  ,  v  , 
col  m  \  cl  y  m’  m’ 


m£ 


X  ,  + 

col 


v  (t-t  )  <  X  <  X*  ft)  implies  the  values  are  P  ,  v  ,  c 
m \  c/  sr'  mramr 
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x  >  X*r(t)  implies  the  values  are  P^,  vf,  p 
c*  Application  as  a  Test  Problem 


(1)  Eulerian  Input 

(a)  Initial  values 


For  X  <  Xg^(0),  the  values  are 


pt>° 


»» » 0 

vt >  0 


For  Xs£(0)  <  X  <  xsr(0)>  the  values  are 


P  >  0 
o 

P  >0 
c 

v  ■  0 
o 


For  X  >  X  (0),  the  values  are 

5 1 


where 


P  >  0 
r 

Pr  >  0 

v  -  0 
r 


P„  >  P  >  P 
%  —  r  o 


P  >  0 
o 

pQ  >  0 

vo  =  ° 


and  one  of  the  quantities  P^,  v^.  p^,  then  all  left  quantities  are  determined. 
Also  give  one  of  the  quantities  Pr,  py,  v  ,  and  all  right  quantities  are 
determined.  Therefore, 


90 


AFWL-TR-67-12  7 


91 


AFWL-TR-67-227 


(3)  Numerical  Values  for  SCTp-VI 
(a)  SCTP-VI-a 

X  =  1  meter 


Xg  (0)  =  75  meters 
Xsr(0)  =  125  meters 

and  let  the  right  hand  boundary  be  at  200  meters  Initially 

?0  =*  104  dynes /cm2 

PQ  -  10-6  gm/cu-3 

v  »  0 
o 

*  108  dynes /cm2 
Pr  *  107  dynes /cm2 

These  values  then  determine  the  remainder  of  the  values  to  be 

P  ^  —  5.997  x  10  8  gm/cm3 
Pj-  *  5.97  x  10  8  gm/cm3 
V£  =  9*13  x  108  cm/sec 
vr  =  -2.88  x  1G8  co/sec 
vs£  ==  1.095  x  107  cm/sec 

Vsr  *  ~3*46  x  1(36  cm/sec 
tc  A  3.468  x  10~4  sec 

Xcol  *  1-13  x  104  cm 

~  3.66  x  !08  dynes/cm2 
Pm£  *  1*43  x  10  5  gm/cm3 
Pmr  -  3.09  x  10- S  gm/cm3 
\  "*  1-96  x  108  cm/sec 
V*Z  =  x  10s  ds/sec 

v*r  =  5.72  x  iO6  cm/sec 

A  reasonable  output  recipe  seems  to  be  to  print  out  at  1  x  10~4  sec,  t  =  3.46 
x  10~4  sec,  and  7  x  10~4  sec.  '  c 
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(b)  SCTP-VI-B 

Same  as  A,  but  P  =  P  =  10s  dynes/cm2.  This  yields  the 
jc  r 

following  values: 

p.  =  p  =  5.997  x  10~6  gta/cm3 
x.  r  i 

=  $.13  x  106  cm/sec» 

vg^  “  -vsr  “  1*095  x  107  cm/sec 

t  *  A  2.282  x  10~4  sec 
col 

X  .  =  l.QQ-x  104  cm 
col 

P  =  7.995  x  10®  dynes/cm2 

ID 

v  =  9.24  x  19~6  cm/sec 
m 

Pm£  *  pmr  “  2,098  x  10~5  gm/cu3 

-v*  =  v*  A  3.65  x  106  cm/scc 
sz  sr 

A  reasonable  output  recipe  seems  to  be  to  print  out  at  1  x  10'4  sec,  t-c  “ 

2.282  x  10-4  sec,  and  7  x  1C-4  sec. 

(4)  Comments  on  the  Computer  Solution 

(a)  General 

The  solution  should  be  correct  until  the  shocks  hit  the  left 
or  right  boundaries  of  the  hydrocode.  Note,  that  SCTP-VI-B  can  be  interpreted 
as  the  reflection  of  a  shock  wave  off  a  wall  where  the  wall  is  at  the  collision 
position.  Because  of  its  symmetry,  SCTP-VI-B  gives  a  code  symmetry  test  on 
shocks . 

(b)  PUFF 

SCTP-VI-A  took  1283  cycles  and  87  seconds  CP  time  on  the 
CDC  6600  computer.  The  major  errors  in  evidence  were  the  spikes  in  the  density 
and  internal  energy.  Hot-thin  spikes  resulted  from  the  initial  discontinuities 
as  we  have  observed  before  in  SCTP-I-A.  A  cold-thick'  spike  resulted  from  the 
shock  collision.  See  figure  6e. 

SCTP-VI-B  took  1500  cycles  and  113  seconds  CDC  6600  CP  time.. 
T- .  major  errors  were  again  in  the  density  and  specific  internal  eneigy.  The 
hot-thin  and  cold-thick  spots  occurred  as  in  A.  No  asymmetries  were  observed. 
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Figure  6f .  SCTP-VI-B  Error  Graph 
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Table 

VI-A 

PUFF  ERROR  TABLE 

FOR  SCTP-VI-A 

Problem  time  =  7  x 

10~4  sec 

PUFF  cycle  =  1283 

Computer  time  =87 

sec 

Number  of  active  zones 

"  201 

Sum  Sqr.  Err. 

Max.  Err. 

Position  of 
Max.  Err. 

Pressure 

0.0392 

+0.352 

X* 

sr 

Velocity 

0.0485 

+0.514 

X* 

sr 

Density 

0.103 

+0.640 

X  n 
col 

Sum  Int.  Energy 

Sum  Kin.  Energy 

Sum  Tot.  Energy 

EXACT 

3.104  x  1012 

1.375  x  1012 

4.480  x  1012 

PUFF 

3.126  x  1012 

1.657  x  1012 

4,782  x  1012 

Table  VI-B 

PUFF  ERROR  TABLE  FOR  SCTP-VI-B 

Problem  time  =  7  x  10" 

-4 

PUFF  cycle  =  1500 

Computer  time  =  113 

Number  of  active  zones 

=  201 

Sum  Sqr.  Err. 

Max.  Err. 

Position  of 
Max.  Err. 

Pressure 

0.0411 

+0.332 

Xsr  Xs\ 

Velocity 

0.0700 

+0.602 

X*  and  X* 
sr  s& 

Density 

0.0907 

-0.516 

and  XDr 

Sum  Int.  Energy 

Sum  Kin.  Energy 

Sum  Tot.  Energy 

EXACT 

7.832  x  1012 

9.430  x  101 1 

8. 775  x  101? 

PUFF 

7.882  x  1012 

8.940  x  10n 

8.776  x  1012 
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7.  trf PROCODE  TEST  PROBLEM  SCXP-VII  ‘ 
a  -  Description  ot'  Problem 

'1V0  shock  waves  are  traveling  in  the  same  direction  which  is  taken  to 
the  right.  When  two  shock  waves  are  traveling  in  the  same  direction,  the  one 
behind  will  always  overtake  the  one  in  front. 

After  overtake  time,  a  rarefaction  travels  back  to  the  left  (for  v  _<  5/3) 
and  a  stronger  shock  travels  on  to  the  right. 

Graphical  representation  is  presented  in  figure  7. 


*  III  1  1 111 

p*’  V 

P£r*  pZr’  V  r 

Pr*  °r ’  vr 

\ 

’s£  v 

sr 

Figure  7a.  Pipe  Plot  for  Initial  Values  in  SCTF-VH 


t 


Figure  7b.  Lagrangian  Wave  Plot  for  SCTP-VII 


■SPSTJ 
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Figure  7c.  3!nifcial  Pressure  Plot  for  SCTP-VII 


A - V 


Figure  7<i.  After  Overtake  Pressure  Plot  for  SCTP-VII 


b.  Derivation  of  Solution 


Aftcr  overtake  v  =>  v  +  ®  (r  )  for  the  shock  traveling  to,  the  right 

,\mr-\nV  (  /  \ 

and  v  =  v£  “  MPJ  for  the  ^rarefaction  traveling  to  the  left.  1$  fP^j  and 

V„(p  f  are  defined  ir.  SCTP-V.j  From  the  above  relations,  v  and  P  are  detei 
&\m}  )  mm 


are  defined  ir.  SCTP- 


From  the  above  relations,  v  and  P  are  deter 


mined*  then  is  determined  from  the  Rankine-Hugoniot  relation  (5:RH) 

,  .  .  P  +P_  , 


~  (?  V  'Pv]  -  — -  (v  -V  ) 
f-1  \  mr  mr  r  r/  2  \  r  mr/ 
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is  determined  from  the  fact  that  the  entropy  does  not  change 
through  a  rarefaction;  therefore. 


Let  v  ,  X  be  the  velocity  and  position  of  the  shock  after  overtake, 
s  s 

Let  v  „,X  „  and  v  ,X  be  the  velocities  and  positions  of  the  left  and  right 
si  si  sr*  sr  *  .  e 

shocks  (defined  only  prior  to  overtake).  Let  (Xo>toJ  be  the  point  where,  over¬ 
take  occurs. 

For  t  <  tQ  and  X  <  XgJi(t),  the  values  are  F^,  v^,  p^ 

For  t  £  tQ  and  X^ft)  <  X  <  Xgr(t),  the  values  are  P£r>  v^,  p£r 

For  t  <  t  and  X  >  X  (t),  the  values  are  P  ,  v  ,  p 
—  o  sr  r  r  r 

Or.  the  other  hand,  for  t  >  tQ,  in  the  region 

x  <  xo +  (-rt)K) 5  V*> 

the  values  are  P^,  v^,  p^.  In  the  region 

XI;U>  <  X  <  Xo  +  ^  vm  +  if  v^(t-to)  -  V t) 

The  velocity  v  goes  linearly  from  v^  on  the  left  up  to  v^  on  the  right  and 


2 


In  the  region 

XD(t)  <  X  <  X  +  v  t-t 
R  o  m  o 
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the  values  are  P  ,  p.  ,  and  v  .  In  the  region 
m  Jim  m 

X  +  y  ft-t  )  <  X  <  X  +  v  (t-t  'l 
o  m  \  o  /  o  s\  of 

the  values  are  p  ,  P  ,  v  .  In  the  region 
mr’  m  m  ’ 

X  >  X  +  v  ft-t  ) 
o  s  \  o/ 

the  values  are  P  ,  p  ,  v  . 

r>  Kr»  r 

e.  Application  as  a  Code  Test  Problem 

(1)  Eulerian  Input 

(a)  Initial  values 

For  X  lxs£(0)»  the  values  are  P£>  v£,  p£ 

For  Xg£(0)  <  X  <_  xgr(0),  the  values  are  P£r>  v£r>  p£r 

For  X  >  X  (0) ,  the  values  are  P  ,  v  ,  p 
sr  ’  r*  r*  r 

These  values  are  determined  when  P  ,  v  ,  p  ,  Pff  ,  and  P  are  given. 

r  r  r  x.ir  * 

(b)  Boundary  values 

At  X  =  0,  hold  the  values  at  P£,  v?,  p£ 

(2)  Lagrangian  Input 

(a)  Initial  values 
Same  as  Eulerian 

(b)  Boundary  values 

At  X(0,t)  =  v£t,  the  values  should  be  P£,  v£,  p£ 

(3)  Numerical  values  for  SCTP-VII 
(a)  SCTP-VII-A 

AX  =  1  meter 

X  .  (0)  =  75  meters 
si 

X  (0)  =  100  meters 
sr 
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pr  =  104  dynes /cm2 
Pr  =  10  6  gra/cm2 


P£r  e  108  dynes/cK2 
Pjj  “  1012  dynes/cm2 
Right  boundary  at  200  meters 

These  values  yield 


c£  =  8*97  x  108  cm/sec 
v£  ~  ^ '•82  x  108  cm/sec 
Pfr  ~  3.60  x  10  8  gm/cm3 


vs  -  4.56  x  108  cm/sec 
vsr  ~  1*095  x  lO'7  cm/sec 
c0  A  1.683  x  10~5  sec 
x0  =  1.018  x  104  an 
Pm  “  ^*32  x  1011  dynes/cm2 
V£r  =  9.13  x  108  cm/sec 
p£r  =  8*997  x  10  6  gm/cm3 
vm  *  5.26  x  10s  cm/sec 
vs  =  6.31  x  108  cm/sec 
pm£  =  1*63  x  10  5  gm/cm3 
pmr  ~  ^*®00  x  10  6  gm/cm3 


A  reasonable  output  recipe  is  print  at  1  x  10"5 
and  3  x  10~s  sec. 


1*653  x  10"5 


(time  of  overtake), 


(b)  SCTP-VII-B 

Same  as  A  except 

AX  =  1  centimeter 
As 1^0)  «  75  centimeters 
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X  (0)  =  100  centimeters 
sr 

Right  boundary  at  200  centimeters 

A  reasonable  output  recipe  is  prirt  at  1  x  10~7,  1.683  x  10“7  (time  cf  overtake) 
and  3  x  10  7  sec. 

(4)  Comments  on  the  Computer  Solution 

(a)  General 

The  computer  solution  should  be  in  agreement  until  the 
rarefaction  wave  or  the  shock  wave  reaches  a  boundary.  The  difference  between 
A  and  B  is  just  the  scaling  of  the  space  mesh.  This  variation  was  introduced 
to  see  what,  if  any,  changes  in  the  computer  solution  that  space  scaling  would 
introduce.  It,  of  course,  should  not  produce  any  changes  in  the  computer  if 
the  code  behaves  properly. 

(b)  PUFF 

There  was  no  discernible  difference  in  the  PUFF  solutions 
of  A  and  B  other  than  the  scale  change.  They  both  took  1642  cycles  and  about 
110  seconds  CP  time  to  run  to  3  x  10“5  sec.  Density- specific  internal  energy 
spikes  are  the  most  noticeable  errors.  Hot-thin  spikes  were  formed  at  the 
initial  shock  positions  and  a  cold^thick  spike  was  formed  at  the  overtake 
position.  See  tables  VII-A  and  VII-B.  In  the  tables,  stands  for  the 
position  of  the  fluid  particle  that  was  at  the  left  hand  shock  front  at  time 
zero. 
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ON  THE  EQUATIONS  MODELING  LINEAR  FLOW  GAS  HYDRODYNAMICS 

We  are  going  to  inspect  the  typical  mathematical  model  of  an  air-like  gas; 
that  is,  a  gas  which  nearly  satisfies  the  following  characteristics: 
homogeneous 

compressible 

inviscid  (nonviscous,  no  internal  friction) 
nonconducting  (no  heat  transfer) 

ideal  or  perfect  equation  of  state  ^PV  =  RT  -  (y-1)  e^ 

The  terms  gas  and  fluid  will  be  used  interchangeably  in  the  following. 

Consider  fluid  flowing  in  a  frictionless,  insulated  pipe.  Let  us  establish 
some  coordinate  systems.  We  will  label  points  in  the  fluid  and  call  this  the 
Lagraugian  coordinate  system.  We  will  label  points  on  the  pipe  and  call  that 
the  Eulerian  coordinate  systam.  Stated  another  way,  the  Lagrangian  coordinate 
system  is  fixed  in  the  fluid  and  the  Eulerian  coordinate  system  is  fixed  on  the 
pipe.  As  to  notation,  x  will  be  used  as  a  label  for  the  fluid  points  and  X 
will  be  used  as  a  label  for  the  pipe  points. 

Let  A  be  the  cross-sectional  area  of  the  pipe,  A  is  taken  as  a  constant 
independent  of  X  and  time.  Also,  che  cross-sectional  shape  is  taken  as  a 
constant  independent  of  X  and  time.  Therefore,  the  pipe  and  the  fluid  may  be 
thought  of  as  one-dimensional  continuums,  i.e.,  topologically  equivalent  to 
some  interval  of  the  real  ii.e,  e.g.,  the  unit  interval  from  0  to  1.  This  is 
because  in  any  cross  section  of  the  fluid  perpendicular  to  the  axis  of  the 
pipe,  all  tha  properties  of  the  fluid  are  constant,  e.g.,  fluid  velocity, 
density,  pressure,  etc.  The  foregoing  defines  what  we  mean  by  linear  flow. 

Consider  the  mass  m  of  fluid  between  the  fluid  points  labeled  x  -  Ax  and 
x  +  Ax.  This  mass  is  constant  because  x  -  Ax  and  x  +  Ax  are  the  labels  of 
fixed  points  in  the  fluid  (not  necessarily  fixed  with  respect  to  the  Dipe)  and 
we  are  demanding  the  conservation  of  mass.  Let  X(x,t)  denote  the  Eulerian 
coordinate  at  time  t  of  the  fluid  point  labeled  x.  Often  X(x,0)  =  x  is  taken 
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as  a  defining  relationship  between  the  fluid  labels  and  the  pipe  labels.  We 
will  follow  this  convention.  Likewise,  p{x,t)  will  denote  the  mass  density  of 
the  fluid  at  time  t  and  at  the  fluid  point  labeled  x. 

p(x,t)  >_  0 

Let 

p(x(x,t),t]  =  p(x,t) 

i.e. ,  p  is  the  density  in  the  Eyl'erian  coordinate  system,  whereas  p  is  the 
density  in  the  Lagrangian  coordinate  system.  The  mass  between  the  fluid  point 
labeled  x  -  Ax  and  x  +  Ax  at  time  t  is  defined  by 

X(x+ex,t) 

m  =  A  ^  p{X,t)dX 
X(x-Ax,t) 

Now  at  time  t  *  0, 

x+Ax 

m  *  A  jj  p<x,0)dx 
x~Ax 


therefore,  by  conservati  a  of  mass, 

x+Ax  X(x+Ax,t) 

^  p(x,0)dx  *  ^  p(X,t)dX 

x-Az  X(x-Asrt) 

Let  us  nmke  a  change  of  variables;  then 


X(x+Ax,t) 

j  p(X,t)dX 
X(x-Ax, t) 


x+Ax 

J  p(s,t)  ~  dx 
x-Ax 
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So,  we  have 


rux  xr ax 

j  p(x„0)dx  =  J  p(x,t)  dx 


Therefore,  under  differentiability  assumptions,  the  integral  equation  can  be 
expressed  by  the  following  differential  equation: 

p(x,0)  =  p(x,t)  (x,t) 

This  is  called  the  conservation  of  mass  egua'  ion. 

Hie  velocity,  v(x,t).  Is  defined  by 


v(x,t)  =  g-jr  (x,t) 


Newton’s  second  law  says  that  the  sum  of  the  external  forces  applied  to  a 
rigid  body  is  equal  to  the  time  rate  of  change  cf  momentum  of  the  rigid  body. 
Consider  the  interval  of  fluid  between  the  fluid  points  labeled  x  +  Ax  and 
x  ~  Ax.  This  fluid  interval  is  not  a  rigid  body.  However,  the  law  still 
applies  to  this  body  if  there  are  no  internal  friccicnai  forces  in  the  fluid. 

In  the  case  of  internal  frictional  forces  in  the  fluid,  we  must  modify  Newton's 
law  to  say  that  the  external  force  applied  to  the  fluid  interval  is  equal  to 
the  internal  frictional  forces  plus  the  time  rate  of  change  of  momentum.  We 
now  postulate  that  our  fluid  model  has  no  internal  frictional  forces.  Such  a 
fluid  is  called  nonviscous  or  inviscid.  The  external  force  (taken  as  positive 
whea  in  the  direction  of  increasing  x)  on  the  interval  at  time  t  is 


A^?(x-Ax,t)  -  P(x+Ax,t)j 


The  momentum  of  the  body  at  time  t  is 


PU,0)  va,t)dt 
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Thus,  Newton’s  second  law  yields 

x+Ax  t+At 

^  p(C»0)  ^v(c,t+At)  -  v(5,t)jdf;  =  {  ^?(x-Ax,t)  -  P(x+Ax,T)jdT 

x-Ax  t 

Under  differentiability  assumptions,  this  integral  equation  may  be  written  as 
the  following  differential  equation: 

P(x,0)  (x,t)  =  -  —  (x,t) 

This  is  called  the  conservation  of  momentum  equation. 

Now  we  will  see  what  relation  the  principle  of  energy  conservation  will 
produce.  Let  E(x,t)  be  the  specific  energy  at  point  x  and  time  t,  i.e., 
x  -  Ax,  x  +  Ax  at  time  t  is 


A 


x+Ax 

J 


x-Ax 


p(x,0)  F.(x,t)dx 


The  conservation  of  energy  principle  states  that  the  increase  in  l;-.e  total 
energy  of  an  interval  is  equal  _o  the  work  dene  on  the  interval  plus  any  energy 
in  the  form  of  heat  which  is  Transferred  to  the  interval.  We  now  postulate 
that  there  will  be  no  heat  transfer  in  our  fluid  model.  Therefore,  energy 
conservation  in  our  model  is  expressed  in  the  following  way: 


x+Ax 

p(x,0)  ^E(x, t+At) 

x-Ax 


t+At 

1 


P(x-Ax,t)  v(x-Ax,t) 


-  P(x+Ax,r)  v(x+Ax,T)dx 

Jnder  differentiability  assumptions,  this  integral  equation  may  be  written  as 
the  following  differential  equation: 
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This  is  called  the  conservation  of  energy  equation.  The  total  energy  is 
composed  of  the  internal  energy  and  the  kinetic  energy.  This  is  expressed  by 

E(x,t)  »  e(x,t)  +  1/2  v2(x,t) 


where  e(x,t)  is  the  specific  internal  energy.  (This  could  be  thought  of  as 
the  defining  equation  for  e.)  This  relation  says  that  the  total  specific 
energy  is  the  sum  of  the  specific  internal  energy  and  the  specific  kinetic 
energy.  Using  the  three  foregoing  conservation  equations,  we  can  derive  the 
following  relation: 


3V 

9t 


The  laws  of  thermodynamics  say  TdS  -  de  +  dW,  and  in  our  fluid  model,  we 
postulate  that  all  work  is  PaV  work  (this  postulate  may  be  provable  from  the 
homogeneous  fluid  postulate).  So  the  equation 


0  = 

9t  9t 


says  that  the  entropy  is  constant  for  each  x.  This  is  called  an  isentropic 
process.  This  process  is  r.ot  necessarily  homentropic.*  We  postulate  the 
so-called  perfect  or  ideal  gas  equation  of  state  PV  =  RT  =  (y-l)e  where  R  and 
Y  are  gas  constants  and  T  is  the  temperature.  If  we  are  not  undergoing  a 
shock  transition  or  passing  through  any  ^ther  discontinuity,  then  PV  =  (y-l)e 
and  TdS  =  de  +  PdV  (the  first  law  of  hydrodynamics)  yield 


*In  general,  S(x,t)  is  constant  in  neither  x  nor  t.  If  it  is  constant  in  t  but 
not  necessarily  in  x,  then  we  say  the  process  is  .isentropic.  If  it  is  constant 
in  both  t  and  x,  .hen  we  say  the  process  is  homenttopic. 
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where 


and 


V  const 


“P 


42 

dT 


P  const 


R  »  (Y-1)CV 


And  if  0  a  da  +  PdV,  then 


V) 

i. 

P 


However,  for  a  fluid  Interval  passing  in  time  At  through  a  shock  traveling  to 
the  right  where  the  pressure  on  the  left  is  P.  and  on  the  right  is  P  ,  we  have 

X.  IT 


t+At 

t 


as  was  shown  in  the  shock  relations  discussion  in  SCTP-I.  Using  the  shock 

relations,  we  can  prove  that  e.  -  e  +  Pflv0  -  P  v  >0.  Therefore,  there  is 

i*  it  /v  x,  r  r 

an  entropy  increase  across  the  shock!  But  we  just  "proved"  by  using  the  three 
conservation  equations  that  the  process  was  isentropic!  How  did  this  contra¬ 
diction  arise?  it  arose  through  the  assumption  that  all  quantities  were  smooth, 
i.e.,  th  the  quantities  in  the  equations 


V  3X 

V  ”  3x 
o 

a2x  M  _  i_  ap 

3t2  p  3x 

o 

35  ,  1  3Pv 

3t  p  3x 

o 


(conservation  of  mass) 


(conservation  of  momentum) 


(conservation  of  energy) 


AFWL-TR-67-127 


were  defined.  B”,t  they  are  not  defined  across  a  .shock  for  there  E,  P,  V,  v 
are  discontinuous  by  the  definition  of  a  shock.  (Notice  that  the  derivation 
did  prove,  however,  that  when  the  flow  is  smooth,  it  is  isentropic.)  Another 
mistake  that  is  sometimes  made  is  in  using  the  first  law  of  thermodynamics, 
dQ  =  de  +  PdV,  and  the  statement  that  the  process  is  adiabatic  (no  heat  transfer) 
to  arrive  at  the  relation  0  =  de  +  PdV.  The  error  in  this  argument  is  that  dO 
is  composed  of  twc  parts,  dQ  internal  +  dQ  external.  Saying  that  the  process 
is  adiabatic  is  saying  that  there  is  no  heat  transfer  to  or  from  the  system, 
i.e.,  dQ  external  =  0.  But  there  may  also  be  internal  friction  in  the  fluid. 

If  there  are  internal  frictional  forces  acting,  then  dQ  internal  >  0,  so  that 
dQ  >  0  when  dQ  external  =  0.  When  dQ  =  0,  the  process  is  isentropic,  for  the 
entropy  is  defined  by  TdS  =  dQ.  But  it  could  happen  that  the  process  is 
adiabatic  but  not  isentropic,  i.e.,  dQ  external  =  0,  dQ  internal  >  0.  And  it 
could  happen  that  the  process  is  isentropic  but  not  adiabatic,  i.e.,  dQ  internal 
>  0,  but  dQ  external  <  0  and  also  dQ  internal  +  dQ  external  =  0.  Thus,  adiabatic 
and  isentropic  are  two  entirely  different  processes.  Our  no-heat-transfer 
postulate  implies  that  we  have  an  adiabatic  process  in  our  model,  i.e.,  dO 
external  «  0.  But  we  have  shown  that  our  model  is  not  isentropic  when  shocks 
occur.  And  we  have  assumed  that  the  fluid  was  in vise id.  The  three  previous 
statements  are  contradictory.  Proof:  The  entropy  increases  across  a  shock; 
therefore,  dQ  >  0.  The  process  is  adiabatic;  therefore,  dQ  external  =  0. 

Hence,  dQ  =  dQ  internal  >  0  across  a  shock.  Therefore,  there  are  internal 
frictional  forces  acting.  This  is  a  contradiction  since  we  assumed  the  fluid 
to  be  inviscidl  Thus,  the  postulates  assumed  for  our  model  are  contradictory. 

I.et  us  review  our  principles  and  postulates: 

Principle  of  conservation  of  mass 

Principle  of  conservation  of  momentum 

Principle  of  conservation  of  energy 

Principle  of  monotone  increasing  entropy 

Postulate  of  homogeneous  fluid 

Postulate  of  comp’  -’usi.ble  fluid 

Postulate  of  nonviscous  fluid 

Postulate  of  nonconducting  fluid 
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Postulate  of  ideal  equation  of  state 
Postulate  of  dW  =  PdV 

To  produce  a  self-consistent  model,  we  must  throw  out  or  modify  some  of  the 
postulates.  The  criteria  for  this  remodeling  will  come  from  a  closer  look  at 
what  we  are  trying  to  model.  So,  let  us  have  a  closer  look  at  air.  The  best 
resolution  of  the  dilemma  seems  to  be  to  recognize  that  the  fluid  we  are  trying 
to  model  has  a  very  small  viscosity  (it  also  has  a  very  small  heat  conductivity 
but  the  more  satisfactory  resolution  seems  to  be  in  the  modification  of  the 
inviscid  postulate),  and  that  this  viscosity  is  not  appreciable  except  at  the 
places  that  32v/3X2  is  large.  The  inclusion  of  viscosity  prevents  the  occurrence 
of  shocks.  The  internal  frictional  force  depends  on  the  product  of  the  vis¬ 
cosity  number  and  82v/3X2.  Thus,  the  entropy  increase  depends  on  this  product. 

As  the  viscosity  number  goes  to  zero,  32v/3X2  can  become  larger  and  larger  and 
apparently  the  limit  of  the  product  is  finite  but  not  zero.  So,  we  effectively 
have  an  ’’implicit  point  viscosity."  That  is,  the  viscosity  effect  (entropy 
increase)  is  felt  only  at  points  where  32v/3X2  i.'  infinite,  i.e.,  at  shocks. 

The  object  of  the  remodeling  is,  therefore,  to  define  the  solution  of  the  limit 
to  be  the  limit  of  the  solution  of  the  equations  including  the  viscosity  as 
the  viscosity  goes  to  zero.  So  we  alter  the  inviscid  postulate  to  "the  fluid 
is  inviscid  except  at  shocks  and  there  the  viscosity  effects  (entropy  increase) 
are  those  described  by  the  shock  relations."  Notice  that  the  erroneous  deriva¬ 
tion  of  0  =  de  +  PdV  can  lead  to  a  success  for  shocks  if  it  is  differenced  like 
the  Rankine-Hugonlot  relation.  That  is,  the  Rankine-Hugoniot  relation  is 


0  «  e. 
% 


e  + 
r 


F.+P 
£.  r 


And  if  the  difference  scheme  used  were 


where  e^ 
be  correc 


ft  n+1 
0  =  eT 


-,n+l 


n  , 
eI  + 


+  P. 


) 


is  the  approximation  to  e(nAt,IAx),  then  the  difference  scheme  will 
t  across  shocks.  But  it  does  not  conserve  the  entropy  in  nonshock 
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regions.*  It  also  does  not  explicitly  conserve  the  total  energy.  Suppose  we 
wish  to  write  the  equations  in  terms  of  the  (X,t)  system  instead  of  the  (x,t) 
system.  That  is,  suppose  the  Eulerian  formulation  is  desired  instead  of  the 
Lagrangian.  The  transformation  from  the  Eulerian  system  to  the  Lagrangian 
system  way  be  expressed  in  this  manner: 


3  ,  (  *  3  .  3X  3  ,  fv  » 

3t  in  3t  +  at  3X  in  ^x,t^ 


1  3  j  /  \  1  3  in  (X,t) 

- -r-  in  (x,t)**— ► - 


p{x,0)  3x 


/3(X,t>  3X 


The  first  relation  follows  from  the  chain  rule  and  expresses  the  so-called 
convection  derivative.  The  second  relation  just  expresses  the  conservation  of 
mass 


The  Lagrangian  formulation  is 


Vo  P 


ax 

3x 


ax 

at 


3v 


1  3P 


3t  p  3x 
0 


3E  B  _  1_  8Pv 

p 

o 


E  =  e  +  1/2  v2 


PV 

Y-l 


%  =  3X 
p  3x 

(conservation  of  mass) 

(definition  of  velocity) 

(momentum  conservation) 

(energy  conservation) 
(definition  of  internal  energy) 
(equation  of  state) 


*As  a  matter  of  fact,  for  pn+^/pn  near  1,  the  entropy  error  is  AS  ~Cy(y-l/2) 
^l-pn+^/pnj.  Thus,  we  have  an  "entropy  increase"  if  pn+*  <  pn  and  an  "entropy 
decrease"  if  pn+^  >  pn. 
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APPENDIX  II 
SOLUTION  METHODS 


1. 

Method  of  Finite  Differences 

JL16 

2. 

Method  of  Lines 

132 

3. 

Method  of  Characteristics 

132 

i 

Analog  Methods 

132 

5.  Hybrid  Methods 
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Method  of  Finite  Differences 

For  illustrative  purposes,  we  are  going  to  use  a  linear  approximation  to 
the  hydrodynamic  equations. 

The  equations  describing  an  isentropic  flow  may  be  written  (see  Appendix  I) 

jo  =  3X 
P  3x 


where 


and 


So 


> 


and 


32X  ^  1_  3? 

3t2  p  3x 

o 

V  =  kpY 


P 


PQ  -  P(x,o) 
S  P(x,o) 


Therefore,  using 


P_  / 3xVy 

Po  “  W 

again,  we  have  « 


Y±1 

jj*  _  _  1_  P_  ^o  /  P_\  Y  32X 

**  po  Po  3x  po  \Po/  3x2 
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Under  the  assumptions  that 


3x 


( y.,o ) 


is  near  zero  and  P/Po  is  near  1,  the  above  equation  is  approximated  by 

3^X  *Po  3^X 

3t2  “  p  3x2 
o 


C2  = 


C  >  0 
o 

is  called  the  ambient  isentropic  sound  speed.  Here 

C  *  C(x,0) 

w 

but  we  will  take  it  to  be  a  constant  for  simplicity  in  the  succeeding  illustra¬ 
tions.  So,  our  equation  for  illustrative  purposes  will  be 


32X 

3^' 


C 


2  32X 

o  3x2 


(2) 


This  linear  second-order  hyperbolic  partial  differential  equation  is  called 
the  wave  equation.  If  central  differences  are  used,  the  difference  equation 
for  the  wave  equation  is 


.n+1 

1 


-  2X,  +  X 


n-1 


+  X 


n 

1-1 


where  Xj  is  to ^approximate  X(IAx,nat).  That  is,  X^  is  the  solution  to  the 
difference  problem  at  the  mash  point  (n,I)  which  is  the  point  (nAt,IAx)  in 
the  tx  plane  and,  therefore,  x"  is  the  difference  solution  which  will  hopefully 
approximate  X(IAx,nAt),  the  differential  solution.  At  this  point,  we  pcse  a 
question  known  as  the  Convergence  Problem:  Will  the  solutions  to  the  differenc 
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problem,  X»,  converge  to  the  solution  to  the  differential  probles,,  xfe.t).  ns 

te*  “  '  The  ans“er  =°  this  proxies,  1»  the  linenr  cnse  depenis  upon  a 
number,  £,  known  as  the  Courant-Friedrichs-Lewy  number 

C  At 

If  i?£  1,  the  answer  is  yes.  Someti»es  jf  is  referred  to  ns  the  CFL  nuwber. 

It  is  called  this  because  the  answer  to  the  Convergence  Problem  in  the  linenr 
case  was  established  in  Courant,  Friedrichs,  and  Lewy's  1928  HATH  A111M.SH  paper. 

John  von  Heunann  made  a  conjecture  that  generalises  this  result  to  the  nonlinear 
equations.  It  is  essentially  that  if 

CAt  , 

AX  —  1 

for  all  x  and  t(c  -  C(x,t)),  then  convergence  would  prevail  in  the  nonlinear 
case  also.  This  conjecture  has  not  been  proved.  Experience  indicates  that 
it  is  reasonable  and  this  is  the  condition  used  in  current  hydrocodes  with 
only  slight  modifications v  The  foregoing  difference  equation  may  be  arranged 
in  the  following  form: 

xf1  .  2*"  -  :<r’-  ♦  p  2(x»+i  -  +  \ 


<3> 

In  this  form,  the  difference  equation  can  be  solved  in  a  "time  marching" 
manner.  That  is,  given  the  values  of  xj,  X*"1  for  L  <  I  <  R,  we  can  solve 
for  the  values  of  xf 1  for  L  +  1  <  1  <  R  -  i.  Suppose  we  "solved"  this 
difference  equation  on  a  machine  which  introduced  round-off  errors.  That  is, 
if  we  have  a  machine  that  has  only  a  finite  number  of  digits  in  its  numbers, 
then  it  will  make  round-off  errors  in  computing 


xi(2~2£2)  -  xr1  ♦  £2(x" 


n  +  xn  \ 
1+1  1-1/ 


t 
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1*1  ^ 

We  denote  this  round-off  error  by  We  denote  the  machine  solution  by  X(n,I). 

We  denote  the  total  deviation  of  X(n,I)  from  X^  by  a”,  i.e., 

Aj  -  X(n,I)  -  Xj 

Now  we  pose  a  question  known  as  the  Stability  Problem:  Will  the  round-off 
errors 


tu  . 

e.  1  <_  m  <_  n 

introduced  in  the  machine  solution  cause  X(n,I)  to  drift  so  far  away  from  X^ 

as  to  be  useless  as  an  approximation?  As  a  first  attack  on  the  stability 

tn 

problem,  we  ask:  If  a  single  perturbation  e^.  is  introduced  at  the  mth  cycle 
and  the  resulting  solution  (with  no  other  errors  or  perturbations  introduced) 
is 


Xj  + 

then  does  Aj(rn)  grow  in  time  (i.e.,  with  increasing  n)?  Sometimes  the  terminol¬ 
ogy  "weak  stability,”  "strong  stability"  is  used  in  regard  to  the  last  two 
questions.  That  is,  if  the  answer  to  the  single  perturbation  question  is  no, 
then  the  difference  equation  solution  is  said  to  be  weakly  stable.  If  the 
answer  to  the  stability  problem  question  is  no,  then  the  difference  equation 
solution  is  said  to  be  strongly  stable.  Let 

A*  £  X(n,I)  -  Xj 

for  n  =  1,  2,  3,  ....  Let  6^^  be  the  propagated  error  due  to  a”  and  A^ 

Then 


x”+1  +  6j+1  .  (x”  +  6”)  (2-2  £2)  -  (x”'1  +  Aj'"1) 


?  \xi+i  +  4m  +  xh  +  4”-i) 


Subtracting  the  unperturbed  equation  (i.e.,  equation  (3)  which  was 
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+  X 


n  \ 
1-1/ 


yields 


If 


are  just  propagated  errors  themselves,  we  could  write 


(4) 


This  is  called  the  equation  of  first  variation  of  equation  (3).  It  is  the 
equation  for  the  propagation  of  perturbations  of  equation  (3).  Notice  that  the 
equation  of  first  variation  of  equation  (3)  is  the  same  equation  as  equation  (3). 
That  is,  the  equation  of  first  variation  of  the  wave  difference  equation  is 
again  the  wave  difference  equation.  This  is  always  true  for  linear  difference 
equations.  Proof:  Let 


L[X+6]  =  L[X]  +  L[6] 

by  linearity.  So  that 

L[X+<5]  -  L[X]  =  L[6] 

Therefore,  the  equation  of  first  variation  is 

L  [  6  ]  =  0 

Note  that  in  the  equation  of  first  variation,  the  introduction  of  r.ew  round-off 

errors  is  not  taken  into  account.  Thus,  including  the  total  round-off  error 
n*f*l 

at  the  last  step,  eT  ,  the  total  error  at  cycle  n+1  is 


.n+1  ,  n+1 

{i  +  ci 


.  n+1 


and  this  is  by  definition  AT  ;  so 

.n+1  ..n+1  ,  n+1 

h  "  si  +  Ei 


That  is,  the  machine  solution  is 
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v  /  ,1  -r\  vn+l  .  xn+l  ,  n+1 
X(n+l,I)  =  Xj  +  °I  +  ei 


Yn+1  ,  .n+1 
X1  +  A1 


Now  we  ask  if  there  is  some  restriction  on  At  and  Ax  such  that  the  propagated 
errors  do  not  grow.  The  answer  to  the  Stability  Problem  question  turns  out  to 
be  Che  same  as  the  Convergence  Problem  question,  i.e.,  ^  1  is  the  condition 

on  At  and  Ax.  That  the  answers  should  be  the  same  is  plausible  from  this  point 
of  view:  The  differential  and  the  machine  solutions  can  both  be  considered  as 
perturbations  of  the  difference  solution.  Then  a  necessary  condition  for  both 
convergence  and  stability  is  that  perturbations  do  not  grow  unboundedly  in  the 
difference  equation.  The  condition  p  <_  1  yields  the  result  that  perturbations 
do  not  grow. 

We  will  now  look  at  John  von  Neumann's  method  of  analyzing  the  perturbation 
growth  or  error  propagation  problem.  Consider  the  propagation  equation 

5»+1  -  + 

The  solutions  to  this  equation  may  be  written  in  terms  of  a  Fourier  series, 

A  typical  term  of  the  series  being 

ikx  at 
a,  e  e 
k 

where  k  is  an  integer.  We  plu|  this  term  into  the  above  difference  equation 
and  solve  for  a  in  terms  of  p  . 

ikAxI  a(n+l)At  _  ikAxI  anAt  (0  9  IkAxI  a(n-l)At 

a,  e  e  -  a,  e  e  o  /  -  a,  e  e 


.  2  (  ikAx(I+l)  anAt  ,  ikAx(I-l)  anAt 

+  Y  i  \  e  e  +  aR  e  e 


eolt  -  (2-2  f)  -  e-°“  +  P2  (eitok  +  <Tikto) 

So  we  have  a  defined  implicitly  as  a  function  of  k  and  y  .  We  can  prove  that 
2  ,2 
if  $  <_  1,  then  the  real  part  of  a  is  5  0.  Therefore,  for  y  ^  1.  the  k 

frequency  component  of  the  solution  to  the  error  propagation  equation  will  not 
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grow  in  time  ^because  the  k  frequency  component  term  is  a^  eikx  eaCj.  This 


holds  for  all  k. 


The  condition  ifi  £  1  implies  ^  ±1  which  is 


C  At 
o 

Ax 


< 


1 


This  condition  is  known  as  the  Courant  or  Courant-Friedrichs-Le.—  or  CFL 
condition.  It  is  the  condition  for  both  convergence  and  stability  of  the  wave 
difference  equation. 

The  condition 


C  At 
o 

Ax 


< 


1 


has  a  geometric  interpretation  in  terms  of  the  domain  of  dependence  of  the 
solutions  to  the  differential  and  difference  equation.  If  the  initial  data 
are  given  at  time  zero  as  X(x,0)  and 


dX  ,  n, 
dt  (X’0) 

then  the  solution  to  the  wave  differential  equation  is  given  by  D'Alembert's 
formula: 


x+C  t 
o 


X(x,t) 


1/2(x(x'co  l>°)  +  X(x+C0  l>0)  +  r  j 


dX 

dt 


(c.o)de 


x-C  t 
o 


Thus,  X(x,t)  depends  on  the  initial  data  given  on  the  interval 

|x-C  t,x+C  t] 

L  0  °  j 

We  may  graphically  describe  the  points  that  X(x,t)  depends  upon  in  figure  A-l, 


1 
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<**■'’  V  K  '■  ,*■ 

o  t  *» 


Domain  oK* 
Dependence 


(x,,t) 


x-C  tjsr.  * 


Figure  A-l.  Domain  of  Dependence  Graph  for  the  Differential  Solution 

Analogously,  the  domain  of  dependence  of  the  difference  solution  x"  is 
graphed  in  figure  A-2  and  is  seen  to  depend  on  the  values 


The  CFL  condition 


requires  that 


XT,  .  for  -  n  <  i  <  n 
I+i  ~  — 


XT,  .  for  -n+l<i<n-l 
I+i  —  — 


C  At 

1 

Ax  — 


>  C 

At  —  o 


This  says  that  the  slope  of  the  bottom  line  should  be  greater  than  or  equal 
to  C  and  the  slope  of  the  top  line  should  be  less  than  or  equal  to  -C  .  Super¬ 
imposing  the  two  graphs  in  the  case  Ax/At  >  CQ  yields  figure  A-3. 
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It  is^seen  that  the  stability  condition  requires  tl  at  the  domain  of  depend¬ 
ence  of  Xj  "overlap"  the  domain  c  dependence  of  X(nAt.IAx).  In  other  words, 

X.£  should  depend  upon  at  least  as  nu.  h  information  (in  the  limit  as  Ax, At  *  0) 
as  X(nAt,IAx). 

At  this  point,  we  will  explore  a  little  farther  into  the  parallel  between 
the  convergence  problem  and  the  stability  problem.  To  reveal  this,  let  us 
consider  x£  as  a  perturbed  solution  of  X(nAt,IAx).  (on  the  other  hand,  X(nAt, 
IAx)  can  be  thought  of  as  a  perturbation  of  xj).  That  is,  in  the  stability 
problem,  the  machine  solution  is  a  perturbed  solution  of  the  finite  difference 
problem  and  the  perturbations  are  the  machine  round-off  errors.  In  the  converg¬ 
ence  problem,  the  perturbations  are  the  deviations  made  in  each  cycle  by  using 
tue  difference  method  to  advance  to  the  next  time  level  instead  of  the 

differential  equation.  To  indicate  the  parallelism,  we  will  redefine  a“  en, 
6j(m)  as  follows:  * 

A“  =  X(nAt.IAx)  -  Xj 

eI  =  the  error  introduced  in  proceeding  to  time  nAt  from  time  (n-l)At 

via  the  difference  equation  rather  than  by  the  differential 
equation 

^(m)  =  the  deviation  due  to  the  propagation  of  e™  by  the  difference 
equation  1 

To  give  a  precise  equation  for  e”  and  we  need  to  select  an  interpolation 

process  in  space  and  time  to  "fill  fa  the  gaps"  between  the  xj.  That  is, 

define  a  function  X*  such  that  x‘  -  ^  (where  t  *  mAt  and  x  -  IAx)  for  n 

and  such  that 


3t 

exists  and  is  an  approximation  to 


At 


end  Xx  is  defined  also  at  all  real  x  between  the  IAx  points. 
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Or»e  of  the  nicest  ways  to  interpolate  in  x  is  to  use  the  difference 
equation  itself  to  fill  in  the  gaps.  That  is,  to  determine  x"At,  set  up  the 
original  vnesh  so  that  for  oome  I,  x  =  lAx,  then  solve  the  difference  equation 
for  Xj,  then  set 


This  can  be  visualized  as  a  translation  of  the  mesh  of  the  difference  solution 
in  the  space  direction.  To  do  something  similar  in  time  could  be  accomplished 
if  we  were  given  a  "band  of  initial  data"  rather  than  a  line  of  initial  data, 
that  is,  initial  data  from  t  =  0  to  some  t  >  0.  Then  we  could  do  a  "time 
direction  translation"  of  the  mesh  to  determine  xljj^.  Perhaps  a  more  practical 
method  (because  we  are  not  normally  given  an  initial  band  of  data)  would  be  to 
vary  the  initial  time  step  in  such  a  way  as  to  produce  a  time  direction  trans¬ 
lation  of  the  mesh.  That  is,  to  define  X^  ,  let  the  initial  time  step  be 

t  -  «  •  GI  (~j 

(where  GI  is  the  greatest  integer  function).  Then  let  At  be  the  timestep 
thereafter. 

£ 

So,  assume  X  is  defined.  Then,  using  D'Alembert's  formula,  the  advance 

,■)  X 

from  Xn  to  Xn  via  the  differential  equation  yields 


x^c  lt  +  x'+clt  + 

\  0  O 


x-CAt 

o  t 

r  sxt 


fc  J 


x-C  At 

O 


e?+1  ■  X*(t+At,x)  -  Xf* 


(where  e  =  nAt  and  x  «  IAx)  and  then  6”+^(n>)  is  the  propagation  of  by 
equation  (4),  We  know  chat  the  condition  for  the  <^(m)  not  to  grow  in  time 
is  that  1.  So  this  is  a  necessary  condition  for  convergence  as  well  as 

for  stability.  In  view  of  the  geometric  interpretation  involving  the  domains 
of  dependence,  this  is  saying  chat  it  is  necessary  that  the  difference  problem 
take  into  account  at  least  as  much  information  (initial  data  domain-wise)  an 
the  differential  problem  if  it  is  to  converge  to  the  differential  solution. 
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A  physical  interpretation  of  the 


C  At 


constraint  is  to  consider  the  difference  equation  as  a  model  of  the  signal 
propagation  process  and  to  interpret  the  constraint  CQAt  £_  Ax  as  due  to  the 
fact  that  the  difference  relation  only  takes  account  of  the  interaction  between 
adjacent  zones  and  so  the  time  step  must  be  restricted  so  that  signals  cannot 
flow  from,  say,  zone  I  to  zone  1+2,  but  that  only  the  signals  from  zones  I  f  1 
and  1  -  1  may  reach  zone  I  in  one  time  step. 

It  turns  out  that  the  condition  ^  £  1  is  also  sufficient  for  convergence. 
For  proof,  see  the  CFL  paper.  For  a  more  general  convergence  proof  for  linear 
problems,  see  Lax's  Equivalence  Theorem  in  Richtmyer's  Difference  Methods  for 
Initial  Value  Problems.  We  will  give  here  the  proof  for  convergence  when  =  1 
because  the  proof  is  illuminating.  Since 

xf 1  ■  xwi  -  + x?-i + (f-1)  (x”+i  - 2x” + x“-i) 

when  p  *  1,  this  becomes 


,n+l 

‘I 


+  X 


n 

1-1 


One  can  see  that  the  solution  is 


where 


2n 


vn  _  Y"  ve  i 
“I  -  jC-j  “l-n+1 
i«0 


(~ 


■j  \ 

i/ 


i 


e  =*  0 
even 


Codd 


-1 


aud  and  X-^  are  the  initial  conditions.  The  foregoing  equation  may  be 
rearranged  to 
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V®  J.  Y®  H-l  v^  „0 

(n  =  lizn  I+n  +  y*  XI-n+2i+2  XI-i 


n+2i  -1 


I-n+2i+l 


By  D'Alembert's  formulc 


x+C  t 


X(t,x)  =  ^x(o,jc-Cot)  +  x(o,x+Cot)j  +  23-  j  (0,x')dx 

°  x-C  t 
o 


This  may  also  be  written 


X(nAt,IAx)  *=  -|(X(0»IAx-nCoAt)  +  x(o,lAx+nCoAt)^ 


IAx+nC  At 
r  0 


JL  (  —  (0 

!Co  J  3t  (0> 


x')dx' 


IAx-nC  At 
o 


Or  using  CoAt„  it  mi.  e  written 


X(nAt.IAx)  -  |jx(of  (I-n)Ax)  +  X  (0,  (Itti) Ax)J 


(I+n) Ax 

+  w-  \  If  ».**>«■ 

o  JK 
(I-n)Ax 


We  also  have 


t  j 

o  .J 


(I+n)  Ax 


(I-n)Ax 
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because  the  2  is  a  X/  AtX,  and  if  we  multiply  by 


we  get 


1  /2to_\ 

i  .x  [Co*tJ 


VC  /  1 


Ax  \  2C 


1  odd 


which  is 


C  At  - 
o  1 

Ax  2C 


o  J  C  At 
o 

£ 


(0,x’)dx’ 


but  in  this  case 


C  At 

t  -TT-l 


So,  as  At, Ax  -*•  0  with  "^  =  1,  then  x"  -*■  X(nAt,lAx). 

Moreover,  the  difference  solution  will  yield  the  exact  solution  to  the  differ¬ 
ential  problem  if  X^  is  defined  properly  when  =  1.  It  is  desired  that 


X°  +  X° 

1+1 _ I-i  x 

2  Aj 


I  x+C  At 
so 


~  (0,x')dx5 

OL 


1  r.-C  nt 
o 


should  approximate 
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'w'' 


If,  3?ter  defining  =  X(0,IAx),  we  define 


IAx+C  At 
o 


-i  =  -1+ilfl-i  _  i_  f  °  8X 
1  ■  2  co.J  .  If  <».*' 


)dx’ 


IAx-C  At 
o 


then 


n-1  xo  o 

V  I-n+2,1+2  AI-n+2j  -1 
Lu  2  I-n+2j+l 

.1=0 


IAx-nC  At 
o 


1_ 

2C 


(  !<»• 


x’)  dx’ 


IAx-C  nAt 
o 


by  the  additivity  of  the  integral.  So  the  difference  solution  is  exactly  the 
same  at  all  points  of  the  grid.  At  this  point,  one  can  make  an  interesting 
observation:  By  the  proofs  in  the  CFL  paper,  the  difference  solution  converges 
to  the  differential  solution  for  t*  <_  1.  But  the  solution  to 


vn+l  vn  vn-l  . 

XI  =  XI+1  "  XI  +  XI-1 


converges  to 


i(x.(°.K  -  |  t 


J  +  x|o. 


\  x  +  7?  t 

,t  +  ft)j  +  ic  |c  If 

'  x  -  4t 

p 


For  <  1,  this  amounts  to  a  "smoothing"  of  the  true  solution.  The  actual 
equation  is 


,n+l  Yn 
I  "  XI+1 


,n-.l 

I 


+  X 


1-1 


+  X 


n 

1-1 
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2  .  ' 

Dropping  the  ^  -  i  A2  XT  term  when  "^  <  1  is  like  adding  its  negative 

which  we  will  denote 

A  A2  x“ 
x  I 


where 


A  = 


=  1  -  -fi1  1  0 


& 

1 

K,  ’ 1 

If  the  solution  to 

k 

,  V 

L- 

?5 

xn+1  =  xn  - 

N; 

I  P-1 

a 

s 

+  A  A2  X 

-  XM  -  X?'1  +  x 1-1  +  2  -  J)  XI+1  -  a?  +  X?-1 


+  A  A2  x“  for  0  <  A  <  1  - ~€  ‘ 
x  I  —  “ 


is  a  continuous  function  of  A*  then  the  solution  varies  from  the  exact  solution 
at  A  =  0  to  an  evidently  more  and  more  smoothed  solution  as  A  grows.  Thus, 
the  A  A2  x”  term  seems  to  have  a  smoothing  effect  reminiscent  of  the  viscosity. 
One  might  guess  that  as  a  function  of  A  of  the  solution  being  converged  to 
would  be  like 


X(t,x)  =  ylxf0,x  - 


h ) + 4 


X  + 


/1=A 


x  + 


/pT 

2C 


C  t 
o 

/I=7 


x  - 


C  t 
o 

/l=A 


||  (0,x*)dx' 


since  "]£  =  /l-A  . 
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2.  Methods  of  Lines  ' 

If,  in  the  differential  equation. 


82X  _  -2  32X 
at*  "  O  8x2 

we  difference  only  in  space,  we  would  have 


32XT (t) 


XI+i(t)  -  2Xx(t)  +  Xi  ;L(t) 


Then  we  have  a  system  of  ordinary  differential  equations.  Getting  approximate 
solutions  to  the  partial  differential  equation  by  solving  this  system  is  called 
the  Method  of  Lines. 

3.  Method  of  Characteristics 

In  Section  II  in  the  discussion  of  ~CTP-II,  the  method  of  characteristics 
is  discussed.  For  a  discussion  of  a  computer  program  to  solve  problems  by 
the  method  of  characteristics,  see  Che  report  SC  47S6(RR),  "'SWAP — A  Computer 
Program  for  Shock  Wave  Analysis."  For  a  discussion  of  the  relative  utility 
of  the  characteristic  method  and  the  fin  diff  method,  see  "On  the  Numerical 
Solution  of  the  Hydro  Equations"  by  Fyfe,  Eng,  and  Toung  (Siam  Review,  Vol  3, 
No.  4,  October  1961),  in  which  the  conclusion  is:  Method  of  characteristics 
is  good  for  relatively  simple  problems.  But  for  problems  not  so  simple,  the 
method  of  finite  differences  is  to  be  preferred.  Moreover,  some  problems  are 
too  complicated  to  handle  by  the  Method  of  Characteristics. 

4.  Analog  Methods 

An  electronic  analog  computer  could  be  used  to  electronically  solve  the 
system  by  the  method  of  lines,  for  example.  Other  analog  techniques  might  be 
based  upon  the  analogy  between  shallow  water  behavior  and  gas  dynamics 
behavior.  The  use  of  wind  tunnels  is  an  Analog  Method  also. 

5  *  Hybrid  Computer  Methods 

A  hybrid  computer  is  a  combination  of  analog  and  digital  computers.  For 
a  report  on  this,  see  AFWL-TR-65-165,  "Development  of  an  Automatic  Device  for 
Solving  Continuum  Mechanics  Problems." 
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APPENDIX  III 

ERROR  FORMATS  AND  ERROR  FUNCTIONS 

The  first  thing  to  be  determined  when  testing  a  code  is  at  precisely  what 
time  and  place  variables  are  computed  by  the  code.  Example:  In  PUFF,  velocity 
is  computed  at  the  zone  boundary,  but  density,  pressure,  and  internal  energy 
are  zone  center  quantities. 

Also,  velocity  is  one  half  of  the  time  step  behind  (timewise)  the  other 
variables  in  PUFF.  When  gradients  are  steep  and  fast  moving,  as  in  shocks, 
errors  of  5  to  10  percent  can  easily  be  made  by  not  plotting  the  variables  at 
the  correct  time  and  place.  Example:  If  one  should  plot  PUFF;s  zone  center 
quantities  at  zone  boundaries  and  plot  PUFF's  velocity  as  if  it  were  at  the 
same  time  as  the  other  variables,  errors  of  5  or  10  percent  may  be  generated. 

It  is  convenient  to  introduce  the  term  'primary  variable.”  By  this,  we 
mean  a  variable  that  the  code  actually  computes  in  its  cycle  advancing  routine. 
Example:  In  PUFF,  the  primary  variables  are  X,  the  zone  boundary  position; 

D,  the  density;  £,  the  stress  or  pressure;  E,  the  internal  energy;  V,  the 
velocity.  Now,  from  these  we  can  compute  other  variables  such  as  momentum, 
kinetic  energy,  total  energy,  etc. 

Those  nonprimary  variables  which  may  be  computed  from  the  primary  variables 
will  be  called  secondary  variables. 

We  wish  to  develop  an  "error  function."  This  error  function  is  to  measure 
the  "distance"  between  the  exact  solution  and  the  code's  solution.  That  is, 
it  is  to  give  us  a  number  or  numbers  indicative  of  the  "closeness"  of  the 
code's  solution  to  the  exact  solution. 

This  question  arises:  What  should  the  error  function  depend  on? 

In  answering  this  question,  the  following  point  should  be  noted:  If  we 
check  the  error,  for  example,  in  momentum  in  PUFF,  we  must  make  an  interpola¬ 
tion  decision  because  momentum  is  the  density-velocity  product  and  PUFF 
computes  these  primary  variables  at  different  places.  For  this  reason,  we 
have  answered  the  error  function  in  this  way:  The  error  function  should 
depend  only  on  primary  variables  at  the  points  in  space  and  time  where  they 


/ 
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are  computed;  that  is,  only  at  zone  boundaries  for  velocities  and  only  at  zone 
centers  for  density,  internal  energy,  and  pressure;  also,  only  at  those  times 
these  variables  are  computed.  This  eliminates  interpolation  decisions. 

Now  actually,  PUFF's  value  at  these  points  is  supposed  to  be  the  average 
over  the  half  zone  to  either  side  of  the  point  (be  it  zone  boundary  or  zone 
center) .  But  the  comparisons  have  all  been  made  with  the  value  of  the  exact 
solution  at  that  point.  Perhaps  it  would  have  been  ..lore  fair  to  PUFF  to  take 
the  exact  solution  and  average  it  over  the  region  composed  of  the  half  zone 
to  the  left  and  the  half  zone  to  the  right  and  then  compare  this  average  with 
PUFF.  But  then  the  difficulty  arises  that  PUFF  may  have  made  an  error  in  its 
zone  boundary  positions.  So,  do  we  average  over  the  region  PUFF  is  averaging 
over  or  over  the  region  the  exact  solution  indicates  we  should  average  over? 

We  have  answered  the  error  function  question  by  providing  several  numbers 
and  graphs.  We  have  provided  the  maximum  error  numbers,  sum  square  error 
numbers,  graphs  indicating  the  characteristic  and  outstanding  distortions, 
etc.  We  normalized  the  maximum  error  number  and  the  sum  square  error  number 
by  dividing  by  the  maximum  absolute  value  of  the  variable  taken  on  in  the 
exact  solution.  The  "point  function"  variables  in  PUFF  which  we  check  are 
the  primary  variables  of  PUFF,  i.e.,  D,  E,  U,  S,  X.  The  other  variables 
checked  were  the  "set  or  interval  functions"  which  were  total  energy,  kinetic 
energy,  and  potential  energy. 

We  now  give  precise  definitions  of  the  sum  square  error  numbers  and  the 
maximum  error  numbers.  Let  fj  be  the  value  of  PUFF  or  whatever  code  is  being 
tested.  Let  fj.  be  the  exact  solutions  value  at  that  place.  Let  N  be  the 
number  of  places  or  zones  which  are  involved. 


sum  square  error 


In  PUFF,  N  is  the  number  of  active  zones  (JSTAR+1  in  FUFF). 


maximum  error 


max 

1<I<N 


max  f  r 

1<I<N  1 


* 


sgn 


'Imax 


Imsx 
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where  I  is  such  that 
max 


I  ~  ,  |  ~  , 

fT  -  fT  =  max  f  -  f 
I  Inax  Imax|  1<I<N  |  I  I| 

(It  is  possible  that  there  could  be  more  than  one  I  .  If  so,  the  program 

chooses  the  larger  one.) 
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APPENDIX  IV 

SUMMARY  OF  PUFF  TEST  RESULTS 

PUFF  s  errors  might  be  categorized  into  overrounds,  underrounds,  overshoots, 
undershoots,  hot-thin  spikes,  and  cold-thick  spikes.  These  may  be  further 
broken  down  into  compression  overrounds,  compression  underrounds,  rarefaction 
overrounds,  and  rarefaction  underrounds.  For  graphical  description  of  these 
distortions,  see  figure  A-r.  The  hot-thin  spikes  occur  where  there  are  shocks 
in  the  initial  flow.  The  cold— thick  spikes  occur  at  shock  collisions  or  shock 
overtakes. 


Figure  A-4.  PUFF  Errors 
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We  will  now  list  some  examples  of  these  errors. 

Compression  overround  in  P,  v,  p,  e  in  SCTP-1,  SCXP-III 
Compression  underround  in  P,  v,  p,  e  in  SCTP-I 
Rarefaction  overround  in  the  pressure  SCTP-V-A 
Rarefaction  tsnderround  in  P,  v,  p,  e  in  SCTP-II 
Rarefaction  overshoot  in  the  velocity  in  SCTP-V-C 
Rarefaction  undershoot  in  P,  v,  p,  e  in  SCTP-II-A 
Hot-thir.  spike  in  SCTP-I,  SCTP-VI,  SCrr-VJI 
Cold-thick  spike  in  SCTP-VI,  SCTP-VII 

The  pressure  and  velocity  are  the  best  behaved  variables  in  PUFF  while 
the  density  and  internal  energy  are  the  worst  behaved. 
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